{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulation-based Calibration in SBI\n",
    "\n",
    "After a density estimator has been trained with simulated data to obtain a posterior, the estimator should be made subject to several diagnostic tests, before being used for inference given the actual observed data. *Posterior Predictive Checks* (see [previous tutorial](https://sbi-dev.github.io/sbi/tutorial/12_diagnostics_posterior_predictive_check/)) provide one way to \"critique\" a trained estimator via its predictive performance. Another important approach to such diagnostics is simulation-based calibration as developed by [Cook et al, 2006](https://www.tandfonline.com/doi/abs/10.1198/106186006X136976) and [Talts et al, 2018](https://arxiv.org/abs/1804.06788).\n",
    "\n",
    "**Simulation-based calibration** (SBC) provides a (qualitative) view and a quantitive measure to check, whether the uncertainties of the posterior are balanced, i.e., neither over-confident nor under-confident. As such, SBC can be viewed as a necessary condition (but not sufficient) for a valid inference algorithm: If SBC checks fail, this tells you that your inference is invalid. If SBC checks pass, this is no guarantee that the posterior estimation is working."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## In a nutshell\n",
    "\n",
    "To run SBC,  \n",
    "\n",
    "1. we sample `theta_o_i` values from the prior of the problem at hand\n",
    "2. we simulate \"observations\" from these parameters: `x_o_i = simulator(theta_o_i)` \n",
    "3. we perform inference given each observation `x_o_i`.\n",
    "\n",
    "This produces a separate posterior $p_i(\\theta | x_{o,i})$ for each of `x_o_i`. The key step for SBC is to generate a set of posterior samples $\\{\\theta\\}_i$ from each posterior (let's call this `theta_i_s`, referring to `s` samples from posterior $p_i(\\theta | x_{o,i})$), and to rank the corresponding `theta_o_i` under this set of samples. A rank is computed by counting how many samples `theta_i_s` fall below their corresponding `theta_o_i` (see section 4.1 in Talts et al.). These ranks are then used to perform the SBC check.\n",
    "\n",
    "### Key ideas behind SBC\n",
    "\n",
    "The core idea behind SBC is two fold: \n",
    "\n",
    "- SBC ranks of ground truth parameters under the inferred posterior samples follow a uniform distribution.  \n",
    "(If the SBC ranks are not uniformly distributed, the posterior is not well calibrated.)\n",
    "\n",
    "- samples from the data averaged posterior (ensemble of randomly chosen posterior samples given multiple distinct observations `x_o`) are distributed according to the prior\n",
    "\n",
    "### What can SBC diagnose?\n",
    "\n",
    "**SBC can inform us whether we are not wrong.** However, it cannot tell us whether we are right, i.e., SBC checks a necessary condition. For example, imagine you run SBC using the prior as a posterior. The ranks would be perfectly uniform. But the inference would be wrong.\n",
    "\n",
    "**The Posterior Predictive Checks (see tutorial 12) can be seen as the complementary sufficient check** for the posterior (only as a methaphor, no theoretical guarantees here). Using the prior as a posterior and then doing predictive checks would clearly show that inference failed. \n",
    "\n",
    "To summarize SBC can:\n",
    "\n",
    "- tell us whether the SBI method applied to the problem at hand produces posteriors that have well-calibrated uncertainties,\n",
    "- and if not, what kind of systematic bias it has: negative or positive bias (shift in the mean of the predictions) or over- or underdispersion (too large or too small variance)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A healthy posterior\n",
    "\n",
    "Let's take the gaussian linear simulator from the previous tutorials and run inference with NPE on it. \n",
    "\n",
    "__Note:__ SBC requires running inference several times. Using SBC with amortized methods like NPE is hence a justified endavour: repeated inference is cheap and SBC can be performed with little runtime penalty. This does not hold for sequential methods or anything relying on MCMC or VI (here, parallelization is your friend, `num_workers>1`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import torch\n",
    "\n",
    "from torch import eye, ones, zeros\n",
    "from torch.distributions import MultivariateNormal\n",
    "\n",
    "from sbi.analysis import check_sbc, run_sbc, get_nltp, sbc_rank_plot\n",
    "from sbi.inference import SNPE, SNPE_C, prepare_for_sbi, simulate_for_sbi\n",
    "from sbi.simulators import linear_gaussian, diagonal_linear_gaussian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "num_dim = 2\n",
    "num_simulations = 10_000\n",
    "\n",
    "prior_mean = ones(num_dim)\n",
    "prior_cov = 2 * eye(num_dim)\n",
    "prior = MultivariateNormal(\n",
    "    loc=prior_mean, covariance_matrix=prior_cov, validate_args=False\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## An ideal case\n",
    "\n",
    "To explore SBC, we make our life easy and assume that we deal with a problem where the likelihood is modelled by an identity mapping and a bit of smear. But to start, we only use an almost vanishing smear of `0.01`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1a68701b2f1144809b8658cfba848f48",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Running 10000 simulations.:   0%|          | 0/10000 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "default_likelihood_loc = 0.0  # let's start with 0 shift\n",
    "default_likelihood_scale = 0.01  # let's smear theta only by a little bit\n",
    "\n",
    "\n",
    "def simulator(theta, loc=default_likelihood_loc, scale=default_likelihood_scale):\n",
    "    \"\"\"linear gaussian inspired by sbibm\n",
    "    https://github.com/sbi-benchmark/sbibm/blob/15f068a08a938383116ffd92b92de50c580810a3/sbibm/tasks/gaussian_linear/task.py#L74\n",
    "    \"\"\"\n",
    "    num_dim = theta.shape[-1]\n",
    "    cov_ = scale * eye(num_dim)  # always positively semi-definite\n",
    "\n",
    "    # using validate_args=False disables sanity checks on `covariance_matrix`\n",
    "    # for the sake of speed\n",
    "    value = MultivariateNormal(\n",
    "        loc=(theta + loc), covariance_matrix=cov_, validate_args=False\n",
    "    ).sample()\n",
    "    return value\n",
    "\n",
    "\n",
    "_ = torch.manual_seed(3)\n",
    "theta, x = simulate_for_sbi(simulator, prior, num_simulations)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "theta: [[1.9352932 1.3774877]]\n",
      "x    : [[1.941461  1.4396194]]\n"
     ]
    }
   ],
   "source": [
    "_ = torch.manual_seed(1)\n",
    "\n",
    "# let's obtain an observation\n",
    "theta_o = prior.sample((1,))\n",
    "x_o = simulator(theta_o)\n",
    "print(\"theta:\", theta_o.numpy())\n",
    "print(\"x    :\", x_o.numpy())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " Neural network successfully converged after 63 epochs."
     ]
    }
   ],
   "source": [
    "_ = torch.manual_seed(2)\n",
    "\n",
    "# we use a mdn model to have a fast turnaround with training.\n",
    "inferer = SNPE(prior, density_estimator=\"mdn\")\n",
    "# append simulations and run training.\n",
    "inferer.append_simulations(theta, x).train();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c09f4fbe3a7c4061921bc431d0f16c91",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Drawing 15000 posterior samples:   0%|          | 0/15000 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "posterior = inferer.build_posterior()\n",
    "posterior_samples = posterior.sample((15_000,), x=x_o)\n",
    "# Generate predictive samples by simulating from posterior samples.\n",
    "posterior_predictive_samples = simulator(posterior_samples)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAHRCAYAAACmZ/R8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAy/UlEQVR4nO3deZCcd33n8c/T19yHRnNIc+gY67J12bKNbQg2joHFOGBzOYEKa7vIZqsi2KVS2ZRXCeuQrIFNKktVvN5NAjZmU2tX1iyCYMDYIUjCgGJAlmRJlqXROTM6ZjSjuXpm+nz2j+6nj0ea++l+np55v6qmMjNPT89XE9yf/v2+v9/vMUzTNAUAgIN8bhcAAFh8CBcAgOMIFwCA4wgXAIDjCBcAgOMIFwCA4wgXAIDjCBcAgOMCbhcALDXv833C7RKABXk1+eKMj2HkAgBwHOECAHAc4QIAcBzhAgBwHOECAHAc4QIAcBxLkd1mmlJsPPV5sFIyDHfrAQAHMHJxW2xc+lJr6sMKGQAocYQLAMBxhAsAwHGECwDAcYQLAMBxhAsAwHGEi8viiaTbJQCA49jn4hLTNPVHLx7Wa8fO6l/dLgYAHEa4uOTZn53V/zvQowrFpXK3qwEAZzEt5oKD3UP6yg/fkiStWV6V+f6R3hG3SgIARxEuRTY8EdNnnz+gWMLUB7eu0O4/eGfm2r/737/UycujLlYHAM4gXIrINE398bcOqefqhFY1VOorH9um8qA/c31oIqbPvfCGixUCgDMIlyL6l+N9+tHRywr5fXr6UztUWx685jHHL41qYCziQnUA4BzCpYj2nx6QJH3s1nZtba+75npnY6r/crB7qJhlAYDjCJciOtQ9LEm6dfWy617f1l4viXABUPoIlyKJJ5J6szcVLjd3XDtqkaRt6dEM4QKg1BEuRXKyb0wTsYSqywLqbKy+7mOskcuh7iElk2YRqwMAZxEuRXIoPRrZ1l4nn+/6d5vc0FKtsoBPI5NxnRkIF7E6AHAW4VIkh3qGJEnbO+qnfEzQ79PWtvTU2PmhwhcFAAVCuBTJwXQzf3t66msqN6fDh74LgFJGuBTBeDSuE+md9zdPM3KRsiMbwgVAKSNciuDohRElkqZaasu0om76Uyqt8Hnr4ogmY4kiVAcAziNcisBq5s80JSZJ7csq1FgdUjxp6uiF4cIWBgAFQrgUgTXFNV0z32IYRmb08gZNfQAlinApAmul2Ez9Fov1uEM9jFwAlCbCpcAGxiLqHpyQpOueJ3Y9N3ekjoc52H21YHUBQCERLgV2OD36uKGp6rqnIF/Pto46GYbUPTjBCckAShLhUmBvzKHfYqktD+qGptQRMSxJBlCKCJcCs1aKzbbfYqHvAqCUES4FZJqmDqeb+dtmsQw516YVNZLEbY8BlCTCpYAGwlFdHY/JMLJhMVvrW9Lh0jdWiNIAoKAIlwI63Z862bitvkLlQf+cfnZ9c6rncvZKWNF40vHaAKCQCJcCOt2fGnV0Nl3//i3TWVlXruqygOJJU2c5fh9AiSFcCuiUFS6NVXP+WcMwtC49ejl5makxAKWFcCkga1rshqa5h4uUunmYpMyJygBQKgiXAjp9xQqXuU+LSdL65lRTv4umPoASQ7gUSDSe1PnBcUnz67lI0vr0yOVkHyMXAKWFcCmQ84PjSiRNVYX8aqktm9dzWMuRz1wJK5ZgxRiA0kG4FIjVzF/bVCXDMOb1HK115aoK+RVLmDrHijEAJYRwKRCrmd/ZOL8pMSm9Yiw9ejnBijEsNYYx9cd0j4UnEC4FYu1xmW8z37Ke5cgAShDhUiDWSrHOeS5DtmygqQ+gBAXcLmCxymygXGC4WMuRGblgScid1jKmee87zeyXEbAdtWQmbV+aU17Lf6A59TXMiJFLAQyGoxoaj0mS1s5jd34ua5f+6StjirNiDECJIFwKwOq3tNVXqDK0sMFh6jlSK8bODow7UR4AFBzhUgCZlWILnBKTJJ8ve8ZYF30XACWCcCmAhRxYeT30XVDSpltSfM2HL/Nh+IzZf/j92Y9gIP8jYPvIuxbM+4BzCJcCOJUZuSxsGbLFOgbmBGeMASgRhEsBnL7izB4XS2Y5MqcjAygRhIvDYomkzg9YB1Y6Oy12uj/MijEAJYF9Lg7rHhxXPGmqIujXitpyR56zrb5CFUG/JmIJnR8cd2y6DXCD4c/ZhzLdXpZrfi7nsf78vSxGMOelLGB7WbO9ITMSicznZjw+5e8zcx537UXbmzz2xFyDkYvDrH7L2sYq+XzOnHPk8xmZUZD1/ADgZYSLwzJnijU7O7rILkemqQ/A+5gWc5jTy5At1uIAwgWeZD+NOGe6y7CP4HOmtAz79FbulJb958qy90UyQiHbc07zPtk+ZRXPme6KRPKvRaLZ3xGL5T9NzhSaOc2MGVIYuTjMmrZaV6iRSz/hAsD7CBcHmaaZGVk4tQzZYoXLqb4xmTQPAXgc4eKgK2NRDU/EZBjOLUO2rF5eKb/P0Fgkrr7RyMw/AAAuoufiIKvf0r6sQuVB/wyPnpuygF+rGip15kpYXX1janFomTPgCNuS4txeSt4yYdn6KqH8I1eMnL6KWZn/v3GzLNtnMe1vi30537CP7O39oGR2GbERyV+K7BvLHg5rhm0rM6M5vyOW/3NmPL8/k39xac40MHJxkDUltq5A+1Bo6gMoFYSLg6yRi9PNfAvLkQGUCqbFHFSoZr4l09RnxRg8LncqLHeqS1LeVJhRkT/1lazP/rcTq7ddC2XfC5u2ZcqJsqk3LPti+dNSvnj2a180f6d98Gq2Np/P9t57LGeaLGnfoT/1NHjeTv8lNEXGyMVBp/oYuQCARLg4JhyJ68LwpKTCjVysFWh9oxGNTE7TQAQAlxEuDrHuPrm8KqRlVaEZHj0/teVBtdSmphgYvQDwMnouDunqT91rxekzxezWNVfr8khEXX1j2rFqWUF/FzCtnCW+1xzxkrv8177cOKfPkqzL3w820Zb9enJZfh8jGcg+Z8L2/i1Wnb2WtF0L2m6DVDac7ZcEx/N7IEYi+3UgYevV5H5hPzE5Z+n1NSctT3e68iLGyMUhp/pSI5dCTYlZrGXONPUBeBnh4pCuAjfzLTfkHAMDAF5FuDjEGknc4PCxL3br2EgJoATQc3FAPJHU2YHCnIZsZz3/+cFxTcYSjh8zA0zJfoxK3jXb8S/lZdf9XMo/xiVRnX9tvDH7khRuy/998Srzup9LUrIue1R+eXU071r4av5+mYrubA+ovN/Wc0nm9Idse1KCub2TaP7vMCYmsz82zd0tlxJGLg44NziuWCJ1a+PWuoqC/q6mmjLVlAeUNJUJNADwGsLFAVb/o7PJuVsbT8UwjMyiAWsRAQB4DdNiDugq8Jliduuaq3Wwe4i+CwprumkwKf9uk/Y7QeaefGy722SyJjtNFWnIXzccqzFyruUfseLvyJ5YvHlFX9619sqhzOchX/601MGB9ryvzxlNmc9Nf/4yaX8kW6svYrsWztbqt90JMzk+kfOF/VTmnL+N/U9q2o+RWTzHwzBycUChT0O2466UALyOcHGAdWvjQm+gtKxZnlqRdp6eCwCPIlwWyDTNgh9YabeqoVJSasUYAHgRPZcFujwS0VgkLp+RuhVxMXQ0pFakXR2PaWQyptry4Aw/ATjAvtw4d/GK/Xj6nOW49uPx4zXZfsVkfX4/Jlqb/dxsyb+d963t3ZnPP9b067xrKwJDmc9jZv7LWqUvf9lwrnNmU97XgbHsf0uBSH5twdHsNf9Q/rXcnpP932uY2a9Nez9mEWPkskBWv2X18iqVBYqz56SmPKiG9OGY3YxeAHgQ4bJAJ/tSp+KtL9KUmKUjPTVGuADwIqbFFuhkeuSyvqW44bKqoVKHuofou8A9eUtsZ38qciL3jpK2wX6sNjtttKJpOO/au5Z1ZT5/d8XFvGvN/ippKC7j7UklRmIya/xKbArJrPMrar6V99jxnGOTLw3V5F2LLsu+JE5O5r/3Lsu5S2XQPkuRs9zasP0tpp0IW0RLj+0IlwXqupwOl+aaGR7prFXpvkv34MQMjwSWgOMTMvamZhFSL/MxBQ5NKnZPpbSqOL1Q5GNabAFM09SJ9LRYsVaKWVgxBqT4h5My9o7KMJX3IVMK7h1X5QhnfbmBcFmAgXBUQ+MxGUbh7+NiR88FSKl4+/rhYU1OtXexH8wNTIstwMn0lFjHskpVhIp7OrE1cum5OqFE0pS/wGeaAdfcbXK6x+Yc/5IM2N7D5vQkkrZXoEQo24NYVp4/5Vvvz76Riub0KoxR2xEqNvXhiJr92dMsVoayvZyayvzlzldD2TeJRsJ+YnLO1/al176c//5tx90okVOfGZu21sWEkcsCdKWnxDYUuZkvSSvrKhTwGYomkro8MjnzDwCLVLxm+tCLVvMy5wb+6gtwMrMzv7jNfEny+wy1L0s19em7YCkLb0gNf+zrrqyvBzawydgNhMsCnLjszh4XSwdNfUDxOp/Me2okQzJzPmRI8XuqFKnlhnpuoOeyAF0u7XGxrKKpD48w7H2G3GuR/Ia7LzZ1j8QXz05xDUzkLyEeSmS/Hjfzp8IiGwMyVtTKfzwqjcaVrPEpujEks86nyWj+y1zuPpfJWP41f84Mc9B26HhwNNsvMaK23kne8S+8Z5cIl3kbDEd1ZSx1ZlGxV4pZWI4MZJl1fsXvqFDMTMz8YBQcETtP1qilrb5CVWXuZDThAsCrGLnMU+ZMMZemxCT2usBl0yxNNnOW39qnkPwT2WkyXzz/jo6BcPY5L/U05F37UeXmzOcx27kxnaHsnSmTtvfMPxtdn/f1rwZWZT4fG8ifeqscz/5+fzR/+s6ITz2dZ+Ye45K0311y+qXSixUjl3k6mTn2xb1wWZU+4v/KWFThCLuQAXgH4TJPmWa+C8uQLbXlQdVXppZZdl9l9ALAOwiXefLCtJiU03cZIFwAeAc9l3kYnojp8kjq2IhiH1hp19FQqcM9wzT14TzbcfD2uygaRvZrM54/LZu3NDmev3rLl9NzCYbz+xFlV7M/Z/rzX56OlLVmPo8m8nsuq6tXZj6fSORvmjx0uS3v65EL2dmG8sv5v6N8IPtvKhvK/zf5wzl3tLQvRY7kHCOTyP/3LqW7T+Zi5DIP1pTYyrpy1bh8i2H2ugDwIsJlHrpcOmb/ejqWsRwZgPcwLTYPJ126Qdj1sNcFBWO/u6R9SW3OLnnTNhWknGkywzZl5pvMTi9V9OdPL5l5y5vz3/smQ2WZz0+YLXnXTocas786Zjst4Gr+cufKS9nnrbiSP2VVMZD9N5YN5p+YbISzpzSbk/mHxZo502TmdHeXXMR3nrRj5DIPJzIHVro/cslMi12dUHKJzu0C8B7CZR66Lrt31L7dyvpy+X2GovGk+kYjM/8AABQB4TJHo5MxXRhODYm9MC0W9PvUWl8uSTo3wB33AHgDPZc5slaKNdeUqa7SG/eJWNdUre7BCb19eVR3dC53uxwsFvb+gK0Hk7vENndZsqS85bhmNJp3yRjP9itC/fnvb42831mWfy2RfWxkrDzvWjLnDpblkfw6Q6P5pVVcyemrXM3vFeWefOwbs92EbzI7M2DaliLn9ZyYnpbEyGXOTrp8zP71bGmrkyQd6R2e4ZEAUByEyxydzNwgzP0pMcvmVitcRlyuBABSCJc58ubIpVZS6s6YkTj3sgDgPnouc2TtcdnQ4p2RS1t9heorgxoaj+nEpTFtba9zuyQsRtfs0cj2Lq65P1fufpUJW+8ih/3Q/twuZrXtjpXlV7NXozX5e1ni5dln8sXz6wyN5RcXGsr2S/yj+bXl3jXTmMxffWnGcvos9n09uY+zX+PIfcxkLBJX71BqI9U6l+4+eT2GYWiLNTV2gb4LAPcRLnNwKj0l1lhdpmVVoRkeXVyb01NjNPUBeAHTYnNwwkObJ+2yIxea+iiSvGmy/KkfM5Zz5It9msh+rEzupZw7WAZs/UP/SPblqiyQ/77Y9E/znLY7SBrj2ekuYzJ/mbRypr7M8Ym8S7knP+f9+6T8qa9rjslZmkuTGbnMQfYGYR4Ml/Ry5LcujiiWWJpzvAC8g3CZA2vkst5DzXzL6oZKVZcFFI0ndap/zO1yACxxhMscnPTwyMXnM3RTq9V3YWoMgLvouczSeDSunqupOVgvjlykVN/l9TODOtI7rI/f2u52OVhK7H2FnLXJpmnrh0SnWdKb28uI5PdDfGU5i2jsvy/3eew9Hb/tCP7c3ontdgC5tSUjtoNgc4+0sR/xktdzWZo9FjtGLrPUlVkpFlKDx1aKWazNlEdZjgzAZYTLLFmbJ71wD5epWE39oxdGuLcLAFcxLTZLJ/qsZcjenBKTpM7GKpUHfRqPJnRmIKwbPLTRE0uYbZrIjOdMi5l+20NzTlq2naZsjtumt3L5Zv8+2ciZNrvmrpHT7MK/Zud93kXezNkxcpmlrsvebeZbAn6fblzJZkoA7iNcZskauazz0GnI12NtpjzKZkoALiJcZiF3pZgXd+fn2sIxMAA8gJ7LLJzuD8s0pYaqkJZXl838Ay6ymvpv9gwrmTTl8019LAbgipz+xDVLepM5y4RtP2b44rlf5F/M+d+5Yb9jpu3Eirzntf/3kZymNvoqc8LIZRaOXUxNMXm532LZ2FKjypBfo5G4utipD8AlhMssvHF+SJJ0c0e9q3XMRsDv07b0/VwOnLvqcjUAlirCZRbeOJ96kb5l1TKXK5mdHek6D5wnXAC4g57LDEYnY3o7fWDljtX17hYzS9lwGXK3EGAmyWn2jlzTO8m9Zj/+Jedx9ueZ7k6Q9t4Nx7g4hpHLDA51D8s0pfZlFWquKXe7nFnZsToVLl19Yxoaj87waABwHuEyA2tqaUeJTIlJqVVtaxurJElvdA+5WwyAJYlwmUE2XOrdLWSObknX+wZNfZQq08z/yLuWzP9IJqb+sD9P7sd0j8WCEC7TSCbNzEoxa6qpVNB3AeAmwmUap6+ENTwRU3kwe2ZXqbDC5WD3kBKckAygyAiXaVhTYtva6hX0l9afauOKGlWF/BqLxHUyfS4aABRLab1iFllmf0uJLEHO5fcZ2p7e9Hng3JCrtQCOmK53As8hXKZhvSiX0kqxXLeuZjMlAHcQLlMYmYxljtkv1XBhpz4AtxAuUzjUPSTTlDoaKtRU4+2TkKdiLUc+3R/W1TCbKQEUD+EyhVKfEpOk+sqQOpuszZSMXgAUD+EyhVLcmX89t6X7LvtOXHG5EgBLCeFyHZOxROa4+lIPl/u3rJQkfe/QBcUS0xzgBwAOIlyu45VjlzUaiautvkKbW0tr86Tdu9c3qrE6pIFwVPtO9LtdDoAlgnC5jhd/1S1J+tit7SV/m+CA36cPb2+TJH37jV6XqwGwVBAuNr1DE3qtK9Wf+MSt7S5X44yP7kiFy6vHLmt4IuZyNQCWAsLF5tu/7pFpSnd2NqijodLtchyxubVWG1qqFY0n9cM3L7pdDoAlgHDJkUyaevHXPZKkT9za4XI1zjEMQx/dkRqFffsAU2MACo9wyfH62UGdHxxXdVlA929d4XY5jnrw5lYZRurf2D047nY5ABY5wiXHi79KjVp+a9tKVYYCLlfjrJV1FXrnDcslSbtp7AMoMMIlbSwS1w/S/YhP3LY4Gvl2H70l9e/a/UavktzjBUABES5p3z3Yq4lYQp1NVSW/cXIqH9iyQlUhv85cCet7hy+4XQ6ARYxwkTQ8HtNXXz0hSfrUO1bJMEp7b8tUqsoC+oN710mSvvLD45qIJlyuCMBiRbhI+ssfHdeVsajWNVfr3961xu1yCuozv7FWbfUVujg8qb/fd9rtcgAsUks+XN44f1XPv35ekvRfH9qiUGBx/0nKg3795w9ukiT97d5Tujg84XJFABajxf1KOoN4Iqk/2X1EppnaxX5n53K3SyqKB7au1G2rl2kiltBfvfy22+UAWISWdLh88xfndOziiOoqgtr1wRvdLqdoDMPQf/nQTZJS5429wZ0qAThsSYZLImnqGz87o7/60XFJ0uP3b1JjdWnebXK+trXX62PpXfuPPfdL/fQkJyYDcM6SC5c3e4b10NM/0xe/d0yTsaTu3dik375t8Rz1Mhd/8sCN2t5Rr6HxmB559nX9rz2nZJrsfwGwcEsmXEYmY/qzfzqqB59+TW/2DqumPKAnP7JFzzxye8kfqz9fDVUh/ePv36nfvq1DSVP6by8f17//h1/rSO8wIQNgQRbXGSfXYZqmdr/Rqy/94LiujEUkSR/e3qo//a0b1VxT7nJ17isP+vWVj23Vto46/dk/HdUrxy7rlWOXtWlFjT5+a7s+fHMrfycAc2aYi+gtqmmamogl1Ht1Qm9dGtXxiyP62akBHeoekiR1NlXpix/erHevb3K30FzRsPSl1tTnuy5IoSrXSnmzZ1h/t++UXjl2WdF46pbIhiHd0lGv929eod/c1KyOZZWqCPldq3ExeJ/vE26XACzIq8kXZ3zMrMPlM8/9csEFLYSpVHik/q8UTyY1EU1oIpbURDSu0cm4hidiil/nzKyKoF//4b71+sxvrPXePhYPhYtleDymfzrUq28d6M0Ec66KoF/Lq0OqCPqVSJqKJ00lkqbKgz5VlQVUFQqoLOiTIWVOO3B74vGZR293uYIswgWlbjbhMutpsR8f71tQMcVUFfJr08pabVpRo00ravS+m1ZoRR1TO7NVVxnUp+9ao0/ftUaXhif16luX9crRS3r9zKAi8aQmYgn1XGXzJYCpzTpc/vJj2wpZx+wYyrwbDvoNVQT9qgj5VR70q7Y8qNqKgOoqgqoI+hft+WDFtqKuXJ++c7U+fedqmaapcDShgbGIBsJRTcYSCvp98hmGfIY0GUtqPBpXOJpQJJZQZgy5aCZeAczWrMPl4duX5nJdZBmGoeqygKrLAlq93P3pOwDe5bEGBABgMSBcAACOI1wAAI4jXAAAjiNcAACOI1wAAI5b9GeLeV6wMrUz3/ocABYBwsVthuGJI18AwElMiwEAHEe4AAAcR7gAABxHuAAAHEe4AAAct6juRAkA8AZGLgAAxxEuAADHES4AAMcRLgAAxxEuAADHzepsMdM0NTo6WuhagIKqqamRYRhulwEsCbMKl9HRUdXV1RW6FqCghoeHVVtb63YZwJIwq30u1xu5jIyMqKOjQ93d3fwHWwD8fZ3HyAUonlmNXAzDmPIFrra2lhe/AuLvC6AU0dAHADiOcAEAOG7e4VJWVqYnnnhCZWVlTtaDNP6+AEoZB1cCABzHtBgAwHGECwDAcYQLAMBxhAsAwHFzDpd9+/bpQx/6kFpbW2UYhr7zne8UoKyl68tf/rJuv/121dTUqLm5WQ899JDefvttt8sCgDmZc7iEw2Ft375dTz/9dCHqWfL27t2rnTt3av/+/Xr11VcVi8X0/ve/X+Fw2O3SAGDWFrQU2TAM7d69Ww899JCDJSFXf3+/mpubtXfvXt19991ulwMAs0LPxeOGh4clSQ0NDS5XAgCzR7h4WDKZ1Oc//3m9613v0pYtW9wuBwBmbVanIsMdO3fu1JEjR/Taa6+5XQoAzAnh4lGf/exn9dJLL2nfvn1qb293uxwAmBPCxWNM09TnPvc57d69W3v27NHatWvdLgkA5mzO4TI2Nqaurq7M12fOnNHBgwfV0NCgVatWOVrcUrRz5049//zz+u53v6uamhpdunRJklRXV6eKigqXqwOA2ZnzUuQ9e/bo3nvvveb7jzzyiJ577jmn6lqyproN7ze+8Q09+uijxS0GAOaJI/cBAI5jKTIAwHGECwDAcYQLAMBxhAsAwHGECwDAcYQLAMBxhAsAwHGECwDAcYQLAMBxhEuRvfDCC6qoqNDFixcz33vssce0bdu2zI3BAKDUcfxLkZmmqZtvvll33323nnrqKT3xxBN69tlntX//frW1tbldHgA4gpFLkRmGoSeffFJf+9rX9OSTT+qpp57Syy+/nAmWl156SRs3btT69ev19a9/3eVqAWB+GLm4ZMeOHTp69KheeeUV3XPPPZKkeDyum266ST/5yU9UV1enW2+9VT//+c+1fPlyl6sFgLlh5OKCl19+WcePH1cikVBLS0vm+6+//ro2b96strY2VVdX6/7779crr7ziYqVACTFNKRpOffCe2XWES5EdOHBADz/8sJ555hndd999+sIXvpC5duHChby+S1tbm3p7e90oEyg9sXHpS62pj9i429UsedzmuIjOnj2rBx54QLt27dInP/lJdXZ26q677tKBAwe0Y8cOt8sDAMcwcimSwcFBfeADH9CDDz6oxx9/XJJ0xx136P7779euXbskSa2trXkjld7eXrW2trpSLwAsBA19D4nH47rxxhu1Z88eGvrAXEXDqSkxSdp1QQpVuVvPEse0mIcEAgH99V//te69914lk0n98R//McECoCQxcgGwODBy8RR6LgAAxxEuABad3W/06m/3nlIyycSMW+i5AFgUkkkz82551+43NaFyrWuq1ntvapn251AYjFwAlLzLI5P6xN//4prvv9F91YVqIBEuABaBr756QscujGS+/k//ZqMk6VA3t7FwC+ECoKT1j0b07QP5xyTdsTa1hP9wz5BYEOsOwgVASfuHX5xVNJHUtvb6zPc2tFSrLODTyGRcZwc4Z8wNhAuAkjURTegf9p+TJD32zjWZ7wf9Pt3UWitJOtQ95EJlIFwAlKxvHejR1fGYOhoq9N4bm/OubU+PZA71DBW/MBAuAEpTImnqmZ+eliR95l1rFfDnv5xt76iTJB3uoanvBsIFQEn657cu6+zAuGrLA/rEbR3XXLd6MEd6hxVLJItcHQgXACXp6+lRy+/euVpVZdfuB1+7vEo1ZQFF4kmduDxa7PKWPMIFQMm5Go7ql2dTGyQfyWnk5/L5DG1jasw1hAuAknM0vWFyzfJKtdSWT/k4a2qMFWPFR7gAKDlv9qZGIpvb6qZ93Pb21PVDjFyKjnABUHKOXEiFxZbWGcKlo16SdOLyqCaiiUKXhRyEC4CScyQ9ctk6w8hlRW25mmrKlEiaOnqB0UsxES4ASsrwREzn0ke6bE7vwp+KYRhMjbmEcAFQUqzTj9vqK7SsKjTj42nqu4NwAVBSZjslZrH6LtYiABQH4QKgpGSa+W3TT4lZrKmzswNhhSPxgtWFfIQLgJJyZJbLkC2N1WVqqimTaUrHL7FTv1gIFwAlYywS1+krYUkzL0POdePK1OjlrYsjMzwSTiFcAJSMty6OyDSzS4xn68aVNZmfR3EQLgBKxps9c+u3WG5i5FJ0hAuAkpFt5s9+SkzKhsvxS6NKJk3H68K1CBcAJeNob2rkMZd+iyStbaxSKODTeDSh84PjhSgNNoQLgJIwEU3oZF9qtddcRy4Bv08bW1J9l2NMjRUF4QKgJLx1aURJM7W0uKV29s18C0394iJcAJQEa3/LlrZaGYYx559nOXJxES4ASoJ1pthMh1VOJRsubKQsBsIFQEk41T8mSdqQ7p3M1Y0rUuHSOzSh4fGYY3Xh+ggXACXhVH9qZ/4NTdXz+vm6yqDa6iskpfo3KCzCBYDnDYajGgxHJUmdTVXzfh6a+sVDuADwvNPpKbHWunJVhgLzfh6a+sVDuADwPKvfckPz/KbELDfR1C8awgWA5y2032KxRi5vXx5VPJFccF2YGuECwPNO9aVHLgvot0jSqoZKVYX8isaTmaP7URiECwDPs4JgoSMXn8/QxhU09YuBcAHgaZF49rDJhfZcJGlTzgnJKBzCBYCnnR8YVyJpqrosoOY53CBsKhvSAXXyMuFSSIQLAE/LrBRrqprXmWJ2G9LTYicujy34uTA1wgWApzm1UsxiHR/TfXVcE9GEI8+JaxEuADwts1LMgX6LlDqyv6EqJNOUuvoYvRQK4QLA06xpsc7GhS1DzrU+HVQn6LsUDOECwLNM08xOizk0cpGyU2Mn+giXQiFcAHhW32hEY5G4fIa0enmlY8+7ocVaMca0WKEQLgA8y5oSW9VQqbKA37HnXW+NXJgWKxjCBYBnOb1SzGJNi/VcnVA4Enf0uZFCuADwLKdXilkaqkJqrA5JYsVYoRAuADwrdwOl09Y3MzVWSIQLAM86XaBpMSmnqc/IpSAIFwCeNB6Nq3doQlJhwoWmfmERLgA86Uz6mP1llUEtqwo5/vxWU5/lyIVBuADwpHMDqWP2Vy93vt8iZafFeocmNMaKMccRLgA8KRsuzm2ezFVfGcoc4c/x+84jXAB40vnB1LTY6obChIuUcwwM4eI4wgWAJ1kjl1UFmhaTpPUt1gGW9F2cRrgA8CTr1saFmhaTGLkUEuECwHOi8aQupJchF3ZajAMsC4VwAeA5vUMTSppSRdCvpnTTvRCsvS6XRiY1Ohkr2O9ZiggXAJ5zbiDVzF/VUCnDMAr2e2rLg2qsToXX2SvjBfs9SxHhAsBzrH7LqgL2WyxrG1O/40w60OAMwgWA52T2uBSw32JZk16NdqafcHES4QLAcwq9gTLXmsZUuJxl5OIowgWA51gbKAu5x8WyNh0u1llmcAbhAsBTTNPM7nEp4rQYIxdnES4APKVvNKLJWFJ+n6G2ZRUF/31r0g39ofGYhsajBf99SwXhAsBTrH5La325gv7Cv0RVhgJaUVsuiakxJxEuADzF2uOyuqHw/RaLNXohXJxDuADwlO4i7nGxWE39s4SLYwgXAJ5yrojNfEtmr8sAu/SdQrgA8JRi7nGxrGHk4jjCBYCnZI5+KWLPJXdazDTNov3exYxwAeAZo5MxDYZTy4GL2XNJHZApjUbiGgizHNkJhAsAz7CmxBqrQ6ouCxTt95YH/WqtS+2pYcWYMwgXAJ6RnRIr3qjFwjEwziJcAHiGNXJxI1ysvS409Z1BuADwjGIeWGnHGWPOIlwAeEYx7+Nil50WY6+LEwgXAJ5h9TusfSfFtIblyI4iXAB4wng0rovDk5KkG5qKHy4dyyrlM6SJWEKXRyJF//2LDeECwBNOp28z3FAVUn1lqOi/PxTwqaOBAyydQrgA8ITT6Rf0ThemxCw09Z1DuADwhNP9Y5KkThemxCycjuwcwgWAJ1jTYjc0VbtWgxUup9JBh/kjXAB4wukr1sjFvXDZ0FIjSTp+adS1GhYLwgWA60zTzIxc3JwW27QiFS49Vyc0Fom7VsdiQLgAcN2lkUmNRxMK+AxXjn6xLKsKqaW2TJL0NqOXBSFcALjOGrWsaqhU0O/uy9LGFbWSpOOXRlyto9QRLgBc54WVYhZraoyRy8IQLgBcdyrTb3GvmW/ZSFPfEYQLANd5YQOlZWPOyIUzxuaPcAHguuy0mPsjl3XN1fL7DA1PxDhjbAEIFwCumowl1Ds0IckbPZfyoD+zmfItmvrzRrgAcNWZK2GZplRbHtDyquIfWHk9G2nqLxjhAsBVmWNfmqtlGIbL1aRsaiFcFopwAeCqTL+l0f1+i8UaubBibP4IFwCuyqwU80C/xXLjytRGylN9Y4olki5XU5oIFwCuskYubtx9cipt9RWqCvkVTSQ5fn+eCBcArsk/sNI702I+n6EN6amxt5gamxfCBYBr+sciGo3E5TOk1cvdO7DyerLHwLAceT4IFwCuOdWXGrW0L6tUWcDvcjX5NrJibEEIFwCu6epLvXB7qd9i2bTSOh2ZcJkPwgWAaw71DEuStrTVuVzJtXJvHDY6GXO5mtJDuABwzeGeIUnStvZ6V+u4nvpKbhy2EIQLAFeEI3F19aWWIW9v997IRZK2tKbqOtI77HIlpYdwAeCKI73DSprSitpyNdeWu13OdW1Nh95hwmXOCBcArjic7rds8+ioRcrW9mYP4TJXhAsAVxxK91u2d9S7Wsd0rIUGXf1jCkfiLldTWggXAK4ohZFLc025VtaVyzSloxfYTDkXhAuAorsajur84LgkaVtbvbvFzGBrevRirWzD7BAuAIrOapCvWV6pusqgy9VML9N3oak/J4QLgKI73D0kyZv7W+y2pmukqT83hAuAojtUAv0WizUtdvpKWCPs1J81wgVA0R0ugZViloaqkNqXVUhiM+VcEC4AiurS8KT6RiPyGdLm1lq3y5kV9rvMHeECoKis/S0bWmpUGQq4W8wsbU2vaKOpP3uEC4Ciyh5W6f1+i4UVY3NHuAAoquzmyXp3C5kD6wDLcwPjGh6nqT8bhAuAokkkzUy4bC+hcKmrDGpN+jbMjF5mh3ABUDQHu4c0PBFTTXlAm1bWuF3OnFj7XQ73DrlaR6kgXAAUzY/fuixJumdDk4L+0nr52dbGirG5KK3/7wIoaf9yvE+S9N4bW1yuZO6spv6vzl2VaZouV+N9hAuAoui5Oq7jl0blM6T3bGxyu5w5u3lVvSqCfvWPRnSc2x7PiHABUBQ/fis1arltdYPqK0MuVzN3ZQG/7uxskCTtO9HvcjXeR7gAKIofp6fE7rux2eVK5u/uDakR176ThMtMCBcABTcWiWv/qQFJ0n0l2G+xWOHyyzNXNR7lzpTTIVwAFNxrJ/sVTSS1ZnmlbmiqcruceetsrFJbfYWiiaT+9fSg2+V4GuECoOCsfstvbmqRYRguVzN/hmFkRi976btMi3ABUFDJpKmfvG0tQS7dfovlng2Nkmjqz4RwAVBQB3uGdGUsqprygG5f2+B2OQv2znWN8vsMnb4SVvfguNvleBbhAqCgfnTkkqTS3JV/PbXlQd2SvskZq8amVvr/nwbgWeFIXC+8fl6S9OHtrS5X45zMkmSmxqZEuAAomH/8ZbdGJuPqbKwqySNfpmKFy8+7BhRLJF2uxpsIFwAFEU8k9cxrZyRJv/fuTvl8pbtKzG5rW53qK4MajcR1sHvI7XI8iXABUBDff/Oieocm1Fgd0kd3tLldjqP8PkP3pEcvX9t32uVqvIlwAeA40zT19+kX3UfuWqPyoN/lipy389518vsMvXLssvakl1oji3AB4LhfnBrQ0Qsjqgj69bt3rna7nILY0FKjR9+5RpL0xe8dUySecLcgjyFcADju79Kjlodva9eyqtI7AXm2/uN716uxukxnroQz/SWkEC4AHPXzrivae6JfPiPVyF/MasuD2vXBTZKkp37cpYvDEy5X5B2ECwDHXA1H9Yf/95Ak6ZPvWKWOhkqXKyq8j9zSpttWL9NELKG/eOkYd6lMI1wAOMI0Te3a/aYujUyqs6lKf/LAjW6XVBSGYeiLD26Wz5B+8OYlfeWHxwkYES4AHPLir3r0wyOXFPQb+pvfuUWVoYDbJRXN5tY6/fmDWySl+k1/9aO3l3zAEC4AFuzMlbD+7HtHJUl/+L6N2tJW53JFxfe7d67WFz+8WZL0P/ec0lf/+aTLFbmLcAGwIL84NaBPfW2/xqMJ3dnZoN+/e3E38afzyDvX6E/T04F/8+OT+ouXji3Z42GWzrgVgKOi8aS++s8n9Ld7T8k0pbWNVfrqb98s/yI65mU+fu/dnUqapr70g+N65rUzOtwzpP/xqR1qqS13u7SiMsylPjEIYNaSSVPHLo5o/+kBfftAr45dHJEk/c7tHfrCb92kqjIX369Gw9KX0icv77oghdy9nfLLRy7qj148rLFIXI3VIf33h2/Wu9c3lvSdOOeCcAEwpf7RiN7sHdLhnmEd7hnWr89d1fBELHO9riKor3x0q+7futLFKtM8Fi6SdLp/TH/wfw7o+KVRSdLyqpDu7FyuO29Yrq1tdbqhqUo15UGXqywMwgUosv2nB9wuIU/SNBVLmIrFk4omkjpzJazDPUN6s2dYF4Ynr3l8Vciv29c26K7O5frILW1q9sp0jwfDRZImogn9+UvH9O0DPYrEr+2/tNSWaW1jlVpqy9VUXaammjJVlwcU9Pnk9xkK+A2F/D6FAqkPv8+Qoezop9gDoTs7l8/qcYQLUGRrHv++2yXMmmFI65qqtbW9Ttva6rS9o15b2uq8eUdJj4aLJRJP6FD3sPafHtDrZwZ14vKo+kYjbpc1Z2e/8sCsHkdDHyiyG5q89aJnGKl3xsGAT0GfoRV15dreXq9t7XXa3Fanajf7KItIWcCvd6xt0DvWNmS+NzIZ0+n+sM4NhNU/GlHfaER9I5MKRxNKJE3FEknFE6n/G00kFY0nFU9mxwNeHhswcgGwOHh85LLUeHBsCwAodYQLAMBxhAsAwHGECwDAcYQLAMBxhAsAwHEsRQawOJimFBtPfR6sLP7WdeQhXAAAjmNaDADgOMIFAOA4wgUA4DjCBQDgOMIFAOA4wgUA4DjCBQDgOMIFAOA4wgUA4DjCBQDgOMIFAOC4gNsFAEuJaZoaHR11uwxgQWpqamTMcDAo4QIU0ejoqOrq6twuA1iQ4eFh1dbWTvsYTkUGimixjFxGRkbU0dGh7u7uGV9kMHde//sycgE8xjAMT75YzFdtbe2i+vd4TSn/fWnoAwAcR7gAABxHuACYs7KyMj3xxBMqKytzu5RFaTH8fWnoAwAcx8gFAOA4wgUA4DjCBQDgOMIFAOA4wgUA4DjCBcCs7du3Tx/60IfU2toqwzD0ne98x+2SFo0vf/nLuv3221VTU6Pm5mY99NBDevvtt90ua94IFwCzFg6HtX37dj399NNul7Lo7N27Vzt37tT+/fv16quvKhaL6f3vf7/C4bDbpc0L+1wAzIthGNq9e7ceeught0tZlPr7+9Xc3Ky9e/fq7rvvdrucOWPkAgAeNDw8LElqaGhwuZL5IVwAwGOSyaQ+//nP613vepe2bNnidjnzwpH7AOAxO3fu1JEjR/Taa6+5Xcq8ES4A4CGf/exn9dJLL2nfvn1qb293u5x5I1wAwANM09TnPvc57d69W3v27NHatWvdLmlBCBcAszY2Nqaurq7M12fOnNHBgwfV0NCgVatWuVhZ6du5c6eef/55ffe731VNTY0uXbokSaqrq1NFRYXL1c0dS5EBzNqePXt07733XvP9Rx55RM8991zxC1pEpron/Te+8Q09+uijxS3GAYQLAMBxLEUGADiOcAEAOI5wAQA4jnABADiOcAEAOI5wAQA4jnABADiOcAEAOI5wAQA4jnABULJeeOEFVVRU6OLFi5nvPfbYY9q2bVvmZltwB8e/AChZpmnq5ptv1t13362nnnpKTzzxhJ599lnt379fbW1tbpe3pHEqMoCSZRiGnnzySX384x/XihUr9NRTT+mnP/1pJlg+8pGPaM+ePbrvvvv0rW99y+VqlxZGLgBK3o4dO3T06FG98soruueeezLf37Nnj0ZHR/XNb36TcCkyei4AStrLL7+s48ePK5FIqKWlJe/ae97zHtXU1LhU2dJGuAAoWQcOHNDDDz+sZ555Rvfdd5++8IUvuF0S0ui5AChJZ8+e1QMPPKBdu3bpk5/8pDo7O3XXXXfpwIED2rFjh9vlLXmMXACUnMHBQX3gAx/Qgw8+qMcff1ySdMcdd+j+++/Xrl27XK4OEiMXACWooaFBx48fv+b73//+912oBtfDajEAi9Z73/teHTp0SOFwWA0NDXrxxRd11113uV3WkkC4AAAcR88FAOA4wgUA4DjCBQDgOMIFAOA4wgUA4DjCBQDgOMIFAOA4wgUA4DjCBQDgOMIFAOA4wgUA4Lj/D95GSCmxAIkaAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 500x500 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# let's do some posterior predictive checks to see if the\n",
    "# posterior predictive samples cluster aournd the observation `x_o`.\n",
    "from sbi.analysis import pairplot\n",
    "\n",
    "fig, ax = pairplot(\n",
    "    samples=posterior_predictive_samples,\n",
    "    points=x_o,\n",
    "    limits=list(zip(x_o.flatten() - 1.0, x_o.flatten() + 1.0)),\n",
    "    offdiag=\"kde\",\n",
    "    diag=\"kde\",\n",
    "    figsize=(5, 5),\n",
    "    labels=[rf\"$x_{d}$\" for d in range(3)],\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The observation `x_o` falls into the support of the predicted posterior samples, i.e. it is within `simulator(posterior_samples)`. Given the simulator, this is indicative that our posterior estimates the data well."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Running SBC\n",
    "\n",
    "We have a working and trained posterior at this point! Hurray! Let's look at the SBC metrics now."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "num_sbc_runs = 1_000  # choose a number of sbc runs, should be ~100s or ideally 1000\n",
    "# generate ground truth parameters and corresponding simulated observations for SBC.\n",
    "thetas = prior.sample((num_sbc_runs,))\n",
    "xs = simulator(thetas)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "SBC is implemented in `sbi` for your use on any `sbi` posterior. To run it, we only need to call `run_sbc` with appropriate parameters. \n",
    "\n",
    "__Note__: For amortized neural posteriors (like in this tutorial), execution of `sbc` is expected to be fast. For posteriors that conduct inference with MCMC and hence are slow, `run_sbc` exposes the use of multiple internal parallel workers to the user. To use this feature, add `num_workers = 2` to the parameters for use of two workers. See the API documentation for details."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "7053a4aca78c4977a1d58f80d271bea5",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Running 1000 sbc samples.:   0%|          | 0/1000 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# run SBC: for each inference we draw 1000 posterior samples.\n",
    "num_posterior_samples = 1_000\n",
    "ranks, dap_samples = run_sbc(\n",
    "    thetas, xs, posterior, num_posterior_samples=num_posterior_samples\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`sbi` establishes two methods to do simulation-based calibration:\n",
    "\n",
    "- metrics to compare the sbc ranks with a uniform distribution \n",
    "- control plots for visual inspections like fig. 1 or 2 in [Talts et al, 2018](https://arxiv.org/abs/1804.06788)\n",
    "\n",
    "The `ranks` count is performed per dimension of `theta`, i.e. on the 1-D marginal posterior estimates. According to theory, the distribution of these ranks (per dimension of `theta`) should turn out to be uniformly distributed. \n",
    "\n",
    "The data average posterior `dap` (see equation 1 of [Talts et al, 2018](https://arxiv.org/abs/1804.06788)) is yet another metric of interest. It is built from singular random samples of the estimated posterior samples for each `xs` above. The `dap` is expected to match the prior distribution used (see equation 1 in [Talts et al, 2018](https://arxiv.org/abs/1804.06788) too). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "check_stats = check_sbc(\n",
    "    ranks, thetas, dap_samples, num_posterior_samples=num_posterior_samples\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `check_stats` variable created contains a dictionary with 3 metrics that help to judge our posterior. The \"first\" two compare the ranks to a uniform distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ranks versus Uniform distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "kolmogorov-smirnov p-values \n",
      "check_stats['ks_pvals'] = [0.28580528 0.12618546]\n"
     ]
    }
   ],
   "source": [
    "print(\n",
    "    f\"kolmogorov-smirnov p-values \\ncheck_stats['ks_pvals'] = {check_stats['ks_pvals'].numpy()}\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Kolmogorov-Smirnov (KS test, see also [here](https://en.wikipedia.org/wiki/Kolmogorov%E2%80%93Smirnov_test#Two-sample_Kolmogorov%E2%80%93Smirnov_test)) as used by `check_sbc` provides p-values `pvals` on the null hypothesis that the samples from `ranks` are drawn from a uniform distribution (in other words `H_0: PDF(ranks) == PDF(uniform)`). We are provided two values as our problem is two-dimensional - one p-value for each dimension. \n",
    "\n",
    "The null hypothesis (of both distributions being equal) is rejected if the p-values fall below a significance threshold (usually `< 0.05`). Therefor, vanishing p-values (`ks_pvals=0`) are interpreted to indicate a vanishing false positive rate to (mistakenly) consider both distrubtions being \"equal\". Both values are above `0.05` and we, therefore, cannot claim that inference clearly went wrong. Nonetheless, we can add additional checks:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "c2st accuracies \n",
      "check_stats['c2st_ranks'] = [0.5765 0.5735]\n"
     ]
    }
   ],
   "source": [
    "print(\n",
    "    f\"c2st accuracies \\ncheck_stats['c2st_ranks'] = {check_stats['c2st_ranks'].numpy()}\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The second tier of metrics comparing `ranks` with a uniform distributions is a `c2st` test (see [here](http://arxiv.org/abs/1610.06545) for details). This is a nonparametric two sample test based on training a classifier to differented one of the ensembles (`ranks` versus samples from a uniform distribution) by being trained on the other. The values reported are the accuracies from cross-validation. If you see values around `0.5`, the classifier was unable to differentiate both ensembles, i.e. `ranks` are very uniform. If the values are high towards `1`, this matches the case where `ranks` is very unlike a uniform distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data averaged posterior (DAP) versus prior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "- c2st accuracies check_stats['c2st_dap'] = [0.474 0.488]\n"
     ]
    }
   ],
   "source": [
    "print(f\"- c2st accuracies check_stats['c2st_dap'] = {check_stats['c2st_dap'].numpy()}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The last metric reported is again based on `c2st` computed per dimension of `theta`. If you see values around `0.5`, the `c2st` classifier was unable to differentiate both ensembles for each dimension of `theta`, i.e. `dap` are much like (if not identical to) the prior. If the values are very high towards `1`, this represents the case where `dap` is very unlike the prior distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visual Inspection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoAAAAHACAYAAAAyfdnSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAlZElEQVR4nO3deZCU1bk/8GeGgWGRYRCURQFREZAgEhdC1GhuUCRq3G5CLG8Uy1yjQBnjlqSMS4wEFzRRyy3JVXK9iYommqiIS9xKVKJEFBdQFMEECW4IiDDMzPn9kZ8dRxZHnJmemfP5VFFlv3367fO+Z/rx26ffpSSllAIAgGyUFrsDAAA0LQEQACAzAiAAQGYEQACAzAiAAACZEQABADIjAAIAZEYABADITFl9GtXW1saSJUuic+fOUVJS0th9AjKUUoqVK1dG7969o7S09X03VUeBxvZZ6mi9AuCSJUuiT58+DdI5gE154403Ytttty12NxqcOgo0lfrU0XoFwM6dOxdWWFFR8fl7BvAJK1asiD59+hTqTWujjgKN7bPU0XoFwI9+rqioqFC4gEbVWn8eVUeBplKfOtr6DrQBAGCTBEAAgMwIgAAAmanXMYDQEtXU1MS6deuK3Q3+vzZt2kRZWVmrPcYPoCURAGmVVq1aFX//+98jpVTsrvAxHTt2jF69ekW7du2K3RWArAmAtDo1NTXx97//PTp27BhbbbWVGadmIKUUVVVV8dZbb8XChQtjwIABrfJizwAthQBIq7Nu3bpIKcVWW20VHTp0KHZ3+P86dOgQbdu2jUWLFkVVVVW0b9++2F0CyJav4LRaZv6aH7N+AM2DagwAkBkBEAAgMwIgNHPjxo2Lww47rNjdAKAVcRII2Xj00EPr3fYrf/pTI/akeTvvvPPijjvuiDlz5myy3QsvvBDnnHNOzJ49OxYtWhS/+MUv4pRTTmmSPgLw+ZgBhEZQVVVV7C40utWrV8f2228fF154YfTs2bPY3QHgMxAAoQHst99+MXHixDjllFOie/fuMXr06IiIuOyyy2Lo0KHRqVOn6NOnT4wfPz5WrVpVeN3UqVOjsrIy7r333hg8eHBsscUWceCBB8abb7650fd66qmnYquttoqLLrpog89XVVXFxIkTo1evXtG+ffvo169fTJ48ufD88uXL47vf/W5stdVWUVFREf/xH/8Rzz77bKE/P/3pT+PZZ5+NkpKSKCkpialTp27wffbYY4+45JJL4tvf/naUl5d/1l0GQBH5CRgayG9/+9s46aSTYubMmYVlpaWlccUVV0T//v3jtddei/Hjx8eZZ54ZV199daHN6tWrY8qUKXHjjTdGaWlp/Nd//Vecfvrp8bvf/W6993jwwQfjiCOOiIsvvjhOOOGEDfbjiiuuiD//+c8xbdq06Nu3b7zxxhvxxhtvFJ7/5je/GR06dIh77rknunTpEtddd1187Wtfi5dffjnGjh0bzz//fMyYMSMeeOCBiIjo0qVLQ+0iGllVVVXU1tYWuxvApygtLS36HZEEQGggAwYMiIsvvrjOso8fE7fddtvFBRdcECeeeGKdALhu3bq49tprY4cddoiIiIkTJ8b555+/3vpvv/32OOaYY+I3v/lNjB07dqP9WLx4cQwYMCD23nvvKCkpiX79+hWee+yxx+Kvf/1rLFu2rDBrN2XKlPjjrbfG/119dRz37W9H26qqKK2p8bNuC1NVVRXz5s2LtWvXFrsrwKcoLy+PQYMGFTUECoDQQHbbbbf1lj3wwAMxefLkmDdvXqxYsSKqq6tjzZo1sXr16ujYsWNE/Ov+uB+Fv4iIXr16xbJly+qsZ9asWXHXXXfFbbfd9qlnBI8bNy7233//GDhwYBx44IFx8MEHxwEHHBAREc8++2ysWrUqunXrVuc1H374YSxcvHhzNptmora2NtauXRtlZWVRVqa0Q3NVXV0da9euLfpsvSoBDaRTp051Hr/++utx8MEHx0knnRSTJk2KLbfcMh577LE4/vjjo6qqqhAA27ZtW+d1JSUlkVKqs2yHHXaIbt26xfXXXx8HHXTQeq/5uC9+8YuxcOHCuOeee+KBBx6Ib33rWzFq1Ki47bbbYtWqVdGrV694+OGH67xm1aJFUdm58+fYepqLsrKyov+0BGxadXV1sbsgAEJjmT17dtTW1sall15auAXatGnTNmtd3bt3jz/+8Y+x3377xbe+9a2YNm3aJkNgRUVFjB07NsaOHRv/+Z//GQceeGC8++678cUvfjGWLl0aZWVlsd122xXar/zYa9u1bRs1jiMDaNWcBQyNZMcdd4x169bFlVdeGa+99lrceOONce211272+rbeeut48MEHY968eXHUUUdt9BvkZZddFjfddFPMmzcvXn755bj11lujZ8+eUVlZGaNGjYqRI0fGYYcdFvfdd1+8/vrr8fjjj8f5l10Wf5s7NyIi+m6zTSz6+99jzpw58fbbb2/0mLKqqqqYM2dOzJkzJ6qqquIf//hHzJkzJxYsWLDZ2whA0zADSDaa+uLOw4YNi8suuywuuuii+PGPfxxf+cpXYvLkyXHMMcds9jp79uwZDz74YOy3335x9NFHx+9///to06ZNnTadO3eOiy++OF555ZVo06ZN7LHHHjF9+vTCLOT06dPjrLPOiuOOOy7eeuut6NmzZ4wcPjy27t49IiIOHT067rzvvvjqV78ay5cvjxtuuCHGjRu3Xl+WLFkSw4cPLzyeMmVKTJkyJfbdd9/1fmIGoHkpSZ882GgDVqxYEV26dIn3338/KioqmqJfsNnWrFkTCxcujP79+0f79u2L3Z0WYeUGZu0677hjg7/PpsamtdeZxt6+NWvWxNy5c6N9+/aOAYRmrKqqKtasWRNDhw5t8P9HfZY64ydgAIDMCIAAAJkRAAEAMiMAAgBkRgCk1arH+U00MWMC0DwIgLQ6H10Wpaqqqsg94ZNWr14dEevf/QSApuU6gLQ6ZWVl0bFjx3jrrbeibdu2hevfsXFVNTXrLVuzZk2DrT+lFKtXr45ly5ZFZWXletcuBKBpCYC0OiUlJdGrV69YuHBhLFq0qNjdaRHWLFu23rL2jXA7uMrKyujZs2eDrxeAz0YApFVq165dDBgwwM/A9fTUJZest2zw1Vc36Hu0bdvWzB9AMyEA0mqVlpa6E0g9pXfeWW+ZfQfQejk4CgAgMwIgAEBmBEAAgMwIgAAAmREAAQAyIwACAGRGAAQAyIwACACQGQEQACAzAiAAQGYEQACAzAiAAACZEQABADJTVuwORERUVVVFbW1tsbsB2aouXf+74Jo1a9ZbVlpaGu3atWuKLgHQiIoeAKuqqmLevHmxdu3aYneFVmbBNdest2zHk04qQk+av3927bresrlz5663rLy8PAYNGiQEArRwRQ+AtbW1sXbt2igrK4uysqJ3h1akbAOzyu3bty9CT5q/+uyr6urqWLt2rdl6gFag2SSusrIyswo0qDYbCCr+xjasvvuqurq6KboDQCNrNgEQgNbnhUmT1ls25KyzitAT4OOcBQwAkBkBEAAgMwIgAEBmBEAAgMwIgAAAmREAAQAyIwACAGRGAAQAyIwACACQGQEQACAzAiAAQGYEQACAzAiAAACZEQABADIjAAIAZEYABADIjAAIAJAZARAAIDMCIABAZgRAAIDMCIAAAJkRAAEAMlNW7A4AwAuTJtV5POSss4rUE8iDGUAAgMwIgAAAmREAAQAyIwACAGRGAAQAyIwACACQGQEQACAzAiAAQGYEQACAzAiAAACZEQABADLjXsAA0Mg+ea/jiKa/33Fj9qG1b19rZAYQACAzAiAAQGYEQACAzDgGcBM293iCxnydYxxaH2MKQFMzAwgAkBkBEAAgMwIgAEBmBEAAgMwIgAAAmREAAQAy4zIwABlpyMsONcdLGDXHPjVnG9pf9WGftnxmAAEAMiMAAgBkRgAEAMiMYwBpkRznAwCbzwwgAEBmBEAAgMwIgAAAmREAAQAyIwACAGRGAAQAyIwACACQGQEQACAzAiAAQGYEQACAzAiAAACZEQABADIjAAIAZEYABADIjAAIAJAZARAAIDMCIABAZgRAAIDMlBW7AwDwSS9MmrTesiFnnVWEnnx2G+o7xdGS/44amxlAAIDMmAGkQXzyW1aO37Ca6zdNYwPAJ5kBBADIjAAIAJAZARAAIDMCIABAZgRAAIDMCIAAAJkRAAEAMiMAAgBkRgAEAMiMO4E0c+4p2bD7oLnerQMaw+Z+dlrS56Q+21jf/dDU2+guPcVhv/+LGUAAgMwIgAAAmREAAQAyIwACAGRGAAQAyIwACACQGQEQACAzAiAAQGYEQACAzAiAAACZcSs4qIeWdGssoGm5ZSctkRlAAIDMCIAAAJkRAAEAMiMAAgBkRgAEAMiMAAgAkBkBEAAgMwIgAEBmBEAAgMwIgAAAmREAAQAy417AALQIjX3P3U+u3/2+P5+WfA/1ltz3+jIDCACQGQEQACAzAiAAQGYEQACAzAiAAACZEQABADIjAAIAZEYABADIjAtB/3+NfYHRYsvhopatfQzroyH3wSfXNeCMMxps3QAUlxlAAIDMmAEEaCUWXHNNlNXWRpva2s/0OrPnkB8zgAAAmREAAQAyIwACAGRGAAQAyIwACACQGQEQACAzAiAAQGZa1HUA63OtqtZ2d4vN5bpezdfmjs0nX+dvHYDNZQYQACAzAiAAQGZa1E/AABSHw0qaj9Z0OEh9/65a8jY2V2YAAQAyIwACAGRGAAQAyIwACACQGQEQACAzAiAAQGYEQACAzLgO4GfUmq6/BADkyQwgAEBmBEAAgMwIgAAAmXEMIABswIbuU+u473/Z3HtDt+R7Sm/uOQCffF1NaWlsN2FCg/Vrc5kBBADIjAAIAJAZARAAIDMCIABAZgRAAIDMOAuYOhrzTict+ewvAGhNzAACAGRGAAQAyIwACACQGQEQACAzTgIBgHpyMhuthRlAAIDMCIAAAJkRAAEAMiMAAgBkRgAEAMiMs4Az5my2f2nM2981JuMHwOYyAwgAkBkBEAAgMwIgAEBmBEAAgMwIgAAAmXEWMAANpj5npzuDvf42d1/Zx/XXkPu4pVxFIsIMIABAdgRAAIDMCIAAAJlpFscALrjmmiirrY02tbWFZZv7O3pz/U2+pd5tYkMc49P47D8AGpMZQACAzAiAAACZEQABADIjAAIAZEYABADIjAAIAJCZZnEZGABg8+Rw2ajmeom3T2pJY2EGEAAgMwIgAEBmBEAAgMxkcQxgczx2oCGPE2hJxxy0BPYnAK2dGUAAgMwIgAAAmREAAQAyIwACAGRGAAQAyIwACACQGQEQACAzAiAAQGYEQACAzDTbO4G0lLsxtJR+0vDqM/b+PgBojswAAgBkRgAEAMiMAAgAkBkBEAAgMwIgAEBmBEAAgMwIgAAAmREAAQAy02wvBA0AsDEutP/5mAEEAMhMtjOAvjkAALkyAwgAkBkBEAAgMwIgAEBmBEAAgMwIgAAAmREAAQAyIwACAGRGAAQAyIwACACQmWzvBEL9uGNK62NMoXny2aQpmQEEAMiMAAgAkBkBEAAgMwIgAEBmBEAAgMwIgAAAmREAAQAyIwACAGRGAAQAyIwACACQGQEQACAzAiAAQGYEQACAzAiAAACZEQABADIjAAIAZEYABADIjAAIAJAZARAAIDMCIABAZgRAAIDMlBW7AwA0nNqSkohS3+2huaotKSl2FyJCAARoFUpLS6NNbW3UlJY2m//BABtWXl4epUX+oiYAArQC7dq1i63efz9SsTsCfKpBgwZFu3btitoHARCglWhTW1vsLgD1UOzwF+EkEACA7AiAAACZEQABADIjAAIAZKbZnATi2lXQvFVXVxe7CwA0kKIHQNeugpahurq6WVy7CoDPr+gB0LWroGUYOnRolJaWNovLFwDw+RQ9AEa4dhW0BO3bty92FwBoIH7LAQDIjAAIAJAZARAAIDMCIABAZgRAAIDMCIAAAJkRAAEAMiMAAgBkRgAEAMiMAAgAkBkBEAAgMwIgAEBmBEAAgMwIgAAAmREAAQAyIwACAGRGAAQAyIwACACQGQEQACAzAiAAQGYEQACAzAiAAACZEQABADIjAAIAZEYABADIjAAIAJAZARAAIDMCIABAZgRAAIDMCIAAAJkRAAEAMiMAAgBkRgAEAMiMAAgAkBkBEAAgMwIgAEBmBEAAgMwIgAAAmREAAQAyIwACAGRGAAQAyIwACACQGQEQACAzAiAAQGYEQACAzAiAAACZEQABADIjAAIAZEYABADIjAAIAJAZARAAIDMCIABAZgRAAIDMCIAAAJkRAAEAMiMAAgBkRgAEAMiMAAgAkBkBEAAgMwIgAEBmBEAAgMwIgAAAmREAAQAyIwACAGRGAAQAyIwACACQGQEQACAzAiAAQGYEQACAzAiAAACZEQABADIjAAIAZEYABADIjAAIAJAZARAAIDMCIABAZgRAAIDMCIAAAJkRAAEAMiMAAgBkRgAEAMiMAAgAkBkBEAAgMwIgAEBmBEAAgMwIgAAAmREAAQAyIwACAGRGAAQAyIwACACQGQEQACAzZfVplFKKiIgVK1Y0Sic+WLeuUdYLNJzG+vx/cv0f1ZvWprHraIRaCi1FY9WBz1JH6xUAV65cGRERffr0+RzdAlq0Ll2a5G1WrlwZXZrovZqSOgoUNHKNq08dLUn1iIm1tbWxZMmS6Ny5c5SUlDRYByP+lVb79OkTb7zxRlRUVDTouike49r6NPaYppRi5cqV0bt37ygtbX1HpzRmHY3wmWuNjGnr1Jjj+lnqaL1mAEtLS2PbbbdtkM5tTEVFhT/wVsi4tj6NOaatcebvI01RRyN85lojY9o6Nda41reOtr6v2QAAbJIACACQmaIHwPLy8jj33HOjvLy82F2hARnX1seYNm/Gp/Uxpq1TcxnXep0EAgBA61H0GUAAAJqWAAgAkBkBEAAgMwIgAEBmih4Ar7rqqthuu+2iffv2MWLEiPjrX/9a7C6xEeedd16UlJTU+Tdo0KDC82vWrIkJEyZEt27dYosttogjjzwy/vnPf9ZZx+LFi+Oggw6Kjh07xtZbbx1nnHFGVFdXN/WmZOvRRx+NQw45JHr37h0lJSVxxx131Hk+pRTnnHNO9OrVKzp06BCjRo2KV155pU6bd999N44++uioqKiIysrKOP7442PVqlV12jz33HOxzz77RPv27aNPnz5x8cUXN/amZU0dbTnU0dahNdTSogbAW265JU499dQ499xz429/+1sMGzYsRo8eHcuWLStmt9iEIUOGxJtvvln499hjjxWe+8EPfhB33nln3HrrrfHII4/EkiVL4ogjjig8X1NTEwcddFBUVVXF448/Hr/97W9j6tSpcc455xRjU7L0wQcfxLBhw+Kqq67a4PMXX3xxXHHFFXHttdfGrFmzolOnTjF69OhYs2ZNoc3RRx8dL7zwQtx///1x1113xaOPPhonnHBC4fkVK1bEAQccEP369YvZs2fHJZdcEuedd1786le/avTty5E62vKooy1fq6ilqYj23HPPNGHChMLjmpqa1Lt37zR58uQi9oqNOffcc9OwYcM2+Nzy5ctT27Zt06233lpY9tJLL6WISE888URKKaXp06en0tLStHTp0kKba665JlVUVKS1a9c2at9ZX0Sk22+/vfC4trY29ezZM11yySWFZcuXL0/l5eXppptuSiml9OKLL6aISE899VShzT333JNKSkrSP/7xj5RSSldffXXq2rVrnTH94Q9/mAYOHNjIW5QndbRlUUdbn5ZaS4s2A1hVVRWzZ8+OUaNGFZaVlpbGqFGj4oknnihWt/gUr7zySvTu3Tu23377OProo2Px4sURETF79uxYt25dnfEcNGhQ9O3btzCeTzzxRAwdOjR69OhRaDN69OhYsWJFvPDCC027Iaxn4cKFsXTp0jpj2KVLlxgxYkSdMaysrIzdd9+90GbUqFFRWloas2bNKrT5yle+Eu3atSu0GT16dMyfPz/ee++9JtqaPKijLZM62rq1lFpatAD49ttvR01NTZ0/4oiIHj16xNKlS4vUKzZlxIgRMXXq1JgxY0Zcc801sXDhwthnn31i5cqVsXTp0mjXrl1UVlbWec3Hx3Pp0qUbHO+PnqO4PhqDTX0mly5dGltvvXWd58vKymLLLbc0zkWgjrY86mjr11JqadnnXgPZGDNmTOG/d9lllxgxYkT069cvpk2bFh06dChizwBaBnWU5qJoM4Ddu3ePNm3arHd20z//+c/o2bNnkXrFZ1FZWRk77bRTLFiwIHr27BlVVVWxfPnyOm0+Pp49e/bc4Hh/9BzF9dEYbOoz2bNnz/VOLqiuro53333XOBeBOtryqaOtT0uppUULgO3atYvddtst/vKXvxSW1dbWxl/+8pcYOXJksbrFZ7Bq1ap49dVXo1evXrHbbrtF27Zt64zn/PnzY/HixYXxHDlyZMydO7fOH/39998fFRUVsfPOOzd5/6mrf//+0bNnzzpjuGLFipg1a1adMVy+fHnMnj270ObBBx+M2traGDFiRKHNo48+GuvWrSu0uf/++2PgwIHRtWvXJtqaPKijLZ862vq0mFraIKeSbKabb745lZeXp6lTp6YXX3wxnXDCCamysrLO2U00H6eddlp6+OGH08KFC9PMmTPTqFGjUvfu3dOyZctSSimdeOKJqW/fvunBBx9MTz/9dBo5cmQaOXJk4fXV1dXpC1/4QjrggAPSnDlz0owZM9JWW22VfvzjHxdrk7KzcuXK9Mwzz6RnnnkmRUS67LLL0jPPPJMWLVqUUkrpwgsvTJWVlelPf/pTeu6559Khhx6a+vfvnz788MPCOg488MA0fPjwNGvWrPTYY4+lAQMGpKOOOqrw/PLly1OPHj3Sd77znfT888+nm2++OXXs2DFdd911Tb69OVBHWxZ1tHVoDbW0qAEwpZSuvPLK1Ldv39SuXbu05557pieffLLYXWIjxo4dm3r16pXatWuXttlmmzR27Ni0YMGCwvMffvhhGj9+fOratWvq2LFjOvzww9Obb75ZZx2vv/56GjNmTOrQoUPq3r17Ou2009K6deuaelOy9dBDD6WIWO/fsccem1L61+ULzj777NSjR49UXl6evva1r6X58+fXWcc777yTjjrqqLTFFlukioqKdNxxx6WVK1fWafPss8+mvffeO5WXl6dtttkmXXjhhU21iVlSR1sOdbR1aA21tCSllD7/PCIAAC1F0W8FBwBA0xIAAQAyIwACAGRGAAQAyIwACACQGQEQACAzAiAAQGYEwAydd955seuuuxa7G5/ZdtttF7/85S8/1zqmTp0alZWVhcctdV8AxdVSa4c6ykcEwBbg4YcfjpKSkvVuEL65Tj/99Dr3KMxZY+2LRx99NA455JDo3bt3lJSUxB133NHg7wHUnzraeBprX0yePDn22GOP6Ny5c2y99dZx2GGHxfz58xv8fXIlAGYkpRTV1dWxxRZbRLdu3T7Xuj5+c+qGaFcsDbEvNuSDDz6IYcOGxVVXXdXg6waKRx1dX2PV0UceeSQmTJgQTz75ZNx///2xbt26OOCAA+KDDz5o8PfKkQDYAPbbb7+YOHFiTJw4Mbp06RLdu3ePs88+Oz5+l7333nsvjjnmmOjatWt07NgxxowZE6+88krh+UWLFsUhhxwSXbt2jU6dOsWQIUNi+vTp8frrr8dXv/rViIjo2rVrlJSUxLhx4yIiora2NiZPnhz9+/ePDh06xLBhw+K2224rrPOjb7z33HNP7LbbblFeXh6PPfbYetP1tbW1cf7558e2224b5eXlseuuu8aMGTMKz7/++utRUlISt9xyS+y7777Rvn37+N3vfrfBfVFSUhLXXHNNfOMb34hOnTrFpEmToqamJo4//vhCPwcOHBiXX355ndeNGzcuDjvssJgyZUr06tUrunXrFhMmTNhk4fvNb34TlZWVm/zmOXXq1Ojbt2907NgxDj/88HjnnXfqPP/JffFRP37+859Hjx49orKyMs4///yorq6OM844I7bccsvYdttt44Ybbtjoe0ZEjBkzJi644II4/PDDN9kO+Bd19N/U0X+ZMWNGjBs3LoYMGRLDhg2LqVOnxuLFi2P27NmbfB311GB3Fc7Yvvvum7bYYov0/e9/P82bNy/93//9X+rYsWP61a9+VWjzjW98Iw0ePDg9+uijac6cOWn06NFpxx13TFVVVSmllA466KC0//77p+eeey69+uqr6c4770yPPPJIqq6uTn/4wx9SRKT58+enN998My1fvjyllNIFF1yQBg0alGbMmJFeffXVdMMNN6Ty8vL08MMPp5T+fbPqXXbZJd13331pwYIF6Z133knnnntuGjZsWKFvl112WaqoqEg33XRTmjdvXjrzzDNT27Zt08svv5xSSmnhwoUpItJ2222X/vCHP6TXXnstLVmyZIP7IiLS1ltvna6//vr06quvpkWLFqWqqqp0zjnnpKeeeiq99tprhf1zyy23FF537LHHpoqKinTiiSeml156Kd15553r7cN+/fqlX/ziFymllC666KLUrVu3NGvWrI2Oy5NPPplKS0vTRRddlObPn58uv/zyVFlZmbp06VJo88l9ceyxx6bOnTunCRMmpHnz5qX/+Z//SRGRRo8enSZNmpRefvnl9LOf/Sy1bds2vfHGG5v4q6i7T26//fZ6tYVcqaP/po5u2CuvvJIiIs2dO7fer2HjBMAGsO+++6bBgwen2trawrIf/vCHafDgwSmllF5++eUUEWnmzJmF599+++3UoUOHNG3atJRSSkOHDk3nnXfeBtf/UQF67733CsvWrFmTOnbsmB5//PE6bY8//vh01FFH1XndHXfcUafNJz+svXv3TpMmTarTZo899kjjx49PKf27cP3yl7/81H0REemUU0751HYTJkxIRx55ZOHxsccem/r165eqq6sLy775zW+msWPHFh5/VLjOPPPM1KtXr/T8889v8j2OOuqo9PWvf73OsrFjx35q4erXr1+qqakpLBs4cGDaZ599Co+rq6tTp06d0k033fSp25mSAAj1oY7+mzq6vpqamnTQQQelvfbaq17t+XRlTTjZ2Kp96UtfipKSksLjkSNHxqWXXho1NTXx0ksvRVlZWYwYMaLwfLdu3WLgwIHx0ksvRUTEySefHCeddFLcd999MWrUqDjyyCNjl1122ej7LViwIFavXh37779/neVVVVUxfPjwOst23333ja5nxYoVsWTJkthrr73qLN9rr73i2Wefrfd6Pq3dVVddFddff30sXrw4Pvzww6iqqlrvrLEhQ4ZEmzZtCo979eoVc+fOrdPm0ksvjQ8++CCefvrp2H777TfZj5deemm9n2BHjhxZ52eZDRkyZEiUlv776IgePXrEF77whcLjNm3aRLdu3WLZsmWbXA/w2aijm26Xcx2dMGFCPP/88/HYY4/Vqz2fzjGAzcR3v/vdeO211+I73/lOzJ07N3bfffe48sorN9p+1apVERFx9913x5w5cwr/XnzxxTrHr0REdOrUqUH6WN/1fLLdzTffHKeffnocf/zxcd9998WcOXPiuOOOi6qqqjrt2rZtW+dxSUlJ1NbW1lm2zz77RE1NTUybNm0ztqB+NtSP+vQNKC51tHXW0YkTJ8Zdd90VDz30UGy77bYN2s+cCYANZNasWXUeP/nkkzFgwIBo06ZNDB48OKqrq+u0eeedd2L+/Pmx8847F5b16dMnTjzxxPjjH/8Yp512Wvz617+OiIh27dpFRERNTU2h7c477xzl5eWxePHi2HHHHev869OnT737XVFREb17946ZM2fWWT5z5sw6ffs8Zs6cGV/+8pdj/PjxMXz48Nhxxx3j1Vdf3ax17bnnnnHPPffEz3/+85gyZcom2w4ePHiD4wI0T+roxuVYR1NKMXHixLj99tvjwQcfjP79+zfJ++bCT8ANZPHixXHqqafG9773vfjb3/4WV155ZVx66aURETFgwIA49NBD47//+7/juuuui86dO8ePfvSj2GabbeLQQw+NiIhTTjklxowZEzvttFO899578dBDD8XgwYMjIqJfv35RUlISd911V3z961+PDh06ROfOneP000+PH/zgB1FbWxt77713vP/++zFz5syoqKiIY489tt59P+OMM+Lcc8+NHXbYIXbddde44YYbYs6cORs9Q+2zGjBgQPzv//5v3HvvvdG/f/+48cYb46mnntrsD/OXv/zlmD59eowZMybKysrilFNO2WC7k08+Ofbaa6+YMmVKHHrooXHvvfd+6s8WDWXVqlWxYMGCwuOFCxfGnDlzYsstt4y+ffs2SR+gpVFHNy7HOjphwoT4/e9/H3/605+ic+fOsXTp0oiI6NKlS3To0KFJ+tCamQFsIMccc0x8+OGHseeee8aECRPi+9//fpxwwgmF52+44YbYbbfd4uCDD46RI0dGSimmT59emBKvqamJCRMmxODBg+PAAw+MnXbaKa6++uqIiNhmm23ipz/9afzoRz+KHj16xMSJEyMi4mc/+1mcffbZMXny5MLr7r777s9cEE4++eQ49dRT47TTTouhQ4fGjBkz4s9//nMMGDCgQfbN9773vTjiiCNi7NixMWLEiHjnnXdi/Pjxn2ude++9d9x9993xk5/8ZKM/8XzpS1+KX//613H55ZfHsGHD4r777ouf/OQnn+t96+vpp5+O4cOHF44jOvXUU2P48OFxzjnnNMn7Q0ukjm5cjnX0mmuuiffffz/222+/6NWrV+HfLbfc0iTv39qVpPSxiyyxWfbbb7/YddddP/ftdQBypY5C0zIDCACQGQEQACAzfgIGAMiMGUAAgMwIgAAAmREAAQAyIwACAGRGAAQAyIwACACQGQEQACAzAiAAQGYEQACAzPw/eslHo8ERd38AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 800x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sbi.analysis import sbc_rank_plot\n",
    "\n",
    "f, ax = sbc_rank_plot(\n",
    "    ranks=ranks,\n",
    "    num_posterior_samples=num_posterior_samples,\n",
    "    plot_type=\"hist\",\n",
    "    num_bins=None,  # by passing None we use a heuristic for the number of bins.\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The two plots visualize the distribution of `ranks` (here depicted in red) in each dimension. Highlighted in grey you see the 99% confidence interval of a uniform distribution given the number of samples provided. In plain english: for a uniform distribution, we would expect 1 out of 100 (red) bars to lie outside the grey area.\n",
    "\n",
    "We also observe, that the entries fluctuate to some degree. This can be considered a hint that `sbc` should be conducted with a lot more samples than 1000. A good rule of thumb is that given the number of bins `B` and the number of SBC samples `N` (chosed to be `1_000` here) should amount to `N / B ~ 20`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqYAAAHACAYAAABwPqpFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABdt0lEQVR4nO3deXydVZ0/8M+z3D252dOkbZqmC22BLlCgFpAWwWF5CcMMo4goiAijLFK2HygqoCL+2JVREAeEEf0pzDgqCApKRdkRZClL9xLatE2z3dz9Wc75/XFz05ZS7nnSPMlN8nm/Xrxsb865z6GWm0/O8j2alFKCiIiIiGiU6aM9ACIiIiIigMGUiIiIiMoEgykRERERlQUGUyIiIiIqCwymRERERFQWGEyJiIiIqCwwmBIRERFRWWAwJSIiIqKyYI72AIry+Tzy+fzg74UQ6OnpQV1dHTRNG8WREREREdEHkVIimUxi8uTJ0PV9n+8sm2B6ww034LrrrhvtYRARERGRR++99x6mTp26z++jlcuVpO+fMU0kEpg2bRree+89xOPxURwZERER0QQgJdbfuBxSCmin/gQzZ88t2aW/vx8tLS3o6+tDVVXVPg+hbGZMQ6EQQqHQHq/H43EGUyIiIiKfCTuPirAOIYFoQ6On/DVc2y7LJpgSERER0fDK5/Po7e1Vamtl+lFcSA9GYn4Oa68YTImIiIjGqWQyifXr1yvNaNrpXjRJCQEd0Uh0BEa3JwZTIiIiovEon0L0me9h1vZNCIfCJZvbdh7dGmAjgIBpjMAA9zSmg6mUEo7jwHXd0R5KWTEMA6ZpsswWERHRRLb5JeDdZ2DmLcAoHTQ1UVjG79ZqRy1DjNlgalkWtm7dikwmM9pDKUvRaBTNzc0IBoOjPRQiIiIaBVu7epDsy2MTJuOPwY8r9ZEBiY5gKw73eWx7MyaDqRACGzduhGEYmDx5MoLBIGcHB0gpYVkWduzYgY0bN2L27NnDUvCWiIiIxpa+RAIGgG69DtviC5X6uEJgcVOAM6ZeWJYFIQRaWloQjY7O5txyFolEEAgE8O6778KyLITDpfeVEBERUfnbunUrksmkUtu+HVtRC4lIOIxvHtOs1Cefz8O2bQbToeBM4N7xz4aIiGicySVgvfkIcr09CIf3rP3+ftHedwAJSDOi/AgpJTRNYzAlIiIiog/x9O1ofPP3qBcShlF6AiqVF+jRNNhGDLZtKz3Ctm0YCgel/MJgSkRERDQWpDohJZCrbIUWqyvZfGvCxmt5ibeDi9CWzSo/JhAIjNrKK4PpPlq+fDkWLVqE22+/HQAwffp0rFixAitWrBjVcREREdH4Ip0cAKCr7RRg6qEl2z+zIYkH+7qxOBrFgQceqPwcwzAYTMeLl156CbHY8F/jdf311+P3v/89Xn31VQSDQfT19Q37M4iIiKiM2VlASkgzDJUdoO5AXdJgwEQkor7PdDQxmA6zhoYGX97Xsix88pOfxNKlS3HPPff48gwiIiIaOVJKZDKZwfvpSwnlUpAApKEWTAdyKQx97JTUHDfBVEqJvCM89wuZuvLJs3Q6jS9/+cv49a9/jcrKSlx++eV7tHn/Ur6mabjrrrvw8MMP48knn0RrayvuvfdeNDQ04Itf/CJeeuklLFy4ED/72c8wc+bMvT77uuuuAwDcd999nv8diYiIqPzkX/klEmtfhOs4Su2r+zoBIaAH1cpAugOB12QwHXl5R+CTdz3nud9DX1qKcEDt9NkVV1yBp556Cr/97W/R2NiIr33ta3jllVewaNGiD+337W9/G7feeituvfVWXHnllfjMZz6DGTNm4Ktf/SqmTZuGL3zhC7jwwgvx2GOPeR4/ERERjUGpTgRe/BHi2az6KXgdgBGEHq2BylScKwBAwlQ4wV8uxk0w9VsqlcI999yDBx54AMcccwwA4P7778fUqVNL9j377LPxqU99CgBw5ZVXYunSpfjGN76B4447DgBw8cUX4+yzz/Zv8ERERFRe8klIKSGMEJKzTwEUV2+teCtEoEKpbXGPKWdMR0HI1PHQl5YOqZ+K9evXw7IsLFmyZPC12tpazJkzp2TfBQsWDP560qRJAID58+fv9loul0N/fz/i8bjq0ImIiGisGjhh75gx9M882ZdHCAlICc6YjgZN05SX5EdaIBAY/HVxP+sHvSaE9z2yRERENAY5eQCA0IPKXd7cnsWf1icHDzWV0pUu7F01xtBtkOMmmPpt5syZCAQCeOGFFzBt2jQAQG9vL9asWYNly5aN8uiIiIhoTHHyAKSnYPrn9Ums7swqn+IHCofD6yvUnzHaGEwVVVRU4JxzzsEVV1yBuro6NDY24uqrrx6xArTt7e3o6elBe3s7XNfFq6++CgCYNWsWKirU9poQERGRfxLd29HX06XUNrxtDeK2A7HLCmopWVtACIET5lShtUbtZH44oOOIeU3KzxhtDKYe3HTTTUilUjjppJNQWVmJyy67DIlEYkSe/c1vfhP333//4O8POuggAMDKlSuxfPnyERkDERER7cXGvyH86NdQk00rT1q5APSgeuH7nFOYKT1q/6n4yNyWoYyy7GnSy3zwCOrv70dVVRUSicQeB4JyuRw2btyItrY2hMNqPzFMNPwzIiIiGkHP/gfyL/8ctm3DNNXOvAip492209A15Vil9rc83YlExsJtn1qARTMn78toh82H5bWh4IwpERER0b5yC3tGd7ScAOvAT5dsbjkC1z25DYnVGrB6m6dHRYPledh7ODCYEhEREe0rx4KUgDDCgF46XvXkbSSsQlWekKleZ7S10kRDZWjIwyx3DKZERERE+8q1AEjAUDvMVCx+XxnSccNxU5T6CCHQ19c3pso/ecVgSkRERPQBLMtSrjFu5jOQUkIqln8aOMfk+VYmTdMG65+PRwymRERERO/jui7WrFmDbDar1H5a51bEXAGYisF0CNeFFs+rM5gSERERTSBu3xaE2p9CCLvf1rg3YScBw9ARjFRAJco6roSEhOvY6OnpURqT4zgIhUIMpkREREQTifGHKzC5cwM0TYOuOqupaZCmWl1SR0hIIRE0DTQ3N8Mw1E7a67qOUIiHn4iIiIgmBimhpQolnHI1c6GpLs+Ha5GrnafU1h3YYxrQNbS0tCgH0/GOwXQfLV++HIsWLcLtt98OAJg+fTpWrFiBFStWjOq4iIiIaIhcG3Lg0NP2RRdBD+974fj3K+4xNXRtxK43HwsYTIfZSy+9hFgsNqzvuWnTJnz729/Gk08+iW3btmHy5Mn47Gc/i6uvvhrBoNpPcURERKTIyQEoHDYSehAqsTFtCTz7bmrw2tBSOpI2AMA0xu9+0aFgMB1mDQ0Nw/6e77zzDoQQ+PGPf4xZs2Zh1apVOPfcc5FOp3HzzTcP+/OIiIjGGyGE8gl7pHsQFC6kZkBTrEv61MYkHl3d73FUEtGAPq4PM3nFYOpBOp3Gl7/8Zfz6179GZWUlLr/88j3avH8pX9M03HXXXXj44Yfx5JNPorW1Fffeey8aGhrwxS9+ES+99BIWLlyIn/3sZ5g5c+YHPvf444/H8ccfP/j7GTNmYPXq1bjzzjsZTImIiEpxLGSeuRt9WzcMLtF/GN3NoSpvwTXCysvsaavwvtOqg2irUVvNdG0LH20b/m0CY9n4CaZSDk69e2KGAcWfVK644go89dRT+O1vf4vGxkZ87WtfwyuvvIJFixZ9aL9vf/vbuPXWW3HrrbfiyiuvxGc+8xnMmDEDX/3qVzFt2jR84QtfwIUXXojHHntMediJRAK1tbXK7YmIiCas9ucQWvX/UJXLwTQVo49hQMTqlWcz7YE9owubIjhuP7WwmUwmEY1yS96uxk8wdXLAvceXbvd+X/gDEChd2iGVSuGee+7BAw88gGOOOQYAcP/992Pq1Kkl+5599tn41Kc+BQC48sorsXTpUnzjG9/AcccdBwC4+OKLcfbZZysPed26dbjjjjs4W0pERKTCShX+J9KIdOtytT4akGk4SPkR9sAxey97RqWUXMZ/n/ETTH22fv16WJaFJUuWDL5WW1uLOXPmlOy7YMGCwV9PmjQJADB//vzdXsvlcujv70c8/uE/ZW3ZsgXHH388PvnJT+Lcc8/1+q9BREQ08QzcY5+LTkZq1im+PMIREkK4yKVT6OmxlfpYllXy+/5EM36CqRkuzH4OpZ/Pdr0xoviT0Qe9Vuo+3o6ODhx99NE4/PDDcffdd/swUiIionHItSElIDT/Yk8hmErUVFVi5swW5X7hsP85ZCwZP8FU05SW5Idq5syZCAQCeOGFFzBt2jQAQG9vL9asWYNly5b59tyiLVu24Oijj8bixYvx05/+lDXPiIiIVLkWAEDq/sUe2y38b7wihsbGRt+eM96Nn2Dqs4qKCpxzzjm44oorUFdXh8bGRlx99dUjEhC3bNmC5cuXo7W1FTfffDN27Ngx+LWmpibfn09ERFRucrkcLMtSahtI9wOuC6mrlX4CgI5+G2u61A9Vd2UcAEDQ5MTRvmAw9eCmm25CKpXCSSedhMrKSlx22WVIJBK+P/eJJ57AunXrsG7duj0OW0mpVsiXiIhoPNm0aRN6e3uV2jZtbkeD60Iz1E/A3/nCDvRmXc/jigR5tei+0GSZJpv+/n5UVVUhkUjssTE4l8th48aNaGtr496MveCfERERjVtS4u2Xn4aVyyAajZZsXr3ht6jY8jf0z/gE+vb7lNIjLnr4PUgJLGiKKJ+0D4kcvnTM/pjcPHFWMz8srw0FZ0yJiIhobHn+R2h94X5IKWEYikvnGpSX8l0hUZy2O2NRLWJBtWf09PQgYHLGdF8wmBIREdGYIre/WfiFbigfaBJmFNn6+aUbArDcnYvJQY932bMu6b5hMCUiIqKxxckDUmLL/Ashpx467G8/GEylgJXLwFYMm67rsmrOPuKfHhEREY0p0i0UsJe6P8vmlishIaHD281M8Xh8tzrl5B1nTImIiGjUSSnVK824FiQAGCGl5n05F3e/2IVkXu2UvSskIIGgqWH27NlKB6yKOGO6b8Z0MC3TggJlgX82REQ0VmSzWaxbtw6O4yi1n9HXA0NKaEYAKt/t1uzIob1PrebprppiJgzDYNgcQWMymBanyTOZDCIR/257GssymQwAcEmBiIjKnm3bSKfTCAaDSkvnurCh6zqMYAQqUba4Z3RWXQj/ekC18phqAg4PM42wMRlMDcNAdXU1Ojs7AQDRaJR/cQZIKZHJZNDZ2Ynq6moYBstWEBFReZNSQgih/P3c0AR0XQMMtckXRxSCaTxsYFq1WpF9ywLyecF8McLGZDAFdl7FWQyntLvq6mpeV0pERGNDahsqEqsRQZVSc00U5klV65LaA8E04GFFXsrCwScG05E1ZoOppmlobm5GY2MjbNse7eGUlUAgwJlSIiIaG3IJRB/+EmZkkjA9FqeXileM2gNL+QEPNUmLwZT7S0fWmA2mRYZhMIQRERGNVekuaMKChAarokW5W75mP4hgpVLbwWCqew+mnDEdWWM+mBIREVF5kVLCddVKM8HKQggJO1iNHUdcr9RFSIntKQduQu2kfV/OhRAC+WwG3d1CqY/ruojH4wymI4zBlIiIiIbV1q1bsW3bNqW2kcQ6TLNtyKD66uevXu/FM++mPY1JSonKWAT77TdbuY9pmlzKH2EMpkRERDSsbNtGMplEdXV1ybaGVlgyN0Pq5R83JwpnS6IBHabivtEAJA6eGkd9fb3yc2jkMZgSERHRsHJdF6FQCOFwuGTbkKlD1zVoiqWfgJ11Sc85pA5zGko/AwD6+vpQX6vWlkYP56eJiIhoWAmhXv9Tk8XST+pzZcVgGvR4yp7L8uWP/w8RERHRsPIUTIX3YGq7hQNMDKbjD5fyiYiIaFi5ruspmEoJCJiDNzSVkh9CXVIADKZjAIMpERERDR8hEOxbj1g6gZAdLdlcT7Rja9LGa8ksfrRts6dHhUxvQZPBtPwxmBIREdFeCSGwefNmOI6j1D624TE0r7ofUkoYRukgaLkSKQFYHg4/AUBThQHdySKdVps1zeVyrEk6BjCYEhER0V45joMdO3YgnU4jFAqVbB/YvhpCCCBSDSsYL9k+Ywu0p1y8HlmKm/5pivK4+nu7EAyoh9m6ujql8dPoYjAlIiKivZJSQkqJmpoaBIOl76aPmoBpGuiZeQqSrR8v2X5jTx63PN2JupCJSEB9qT1jGGhpaUFNTY1yHyp/3GxBREREeyWEGLw3XoXuZgv9TLWaocUDTx63iwIAl+bHIc6YEhER0V5pG1Zi8upHEAgElQ4PBfvfBQBIQy2YDhywh6l7K/0EMJiORwymREREtFeBZ29DVbIXhmHASw50IrVK7dyBGVPDYzDVNI3BdBxiMCUiIqIP5tqAkwMA9M78V2im2uEhJ1IHKz5Dre1gMPU2NE3TWP5pHGIwJSIimkBs20Z/f79SW81KIjZQJirRdiJ0s/ThJ6+cwiVOMDzMfnLGdPxiMCUiIppA+vr6sG7dOqW2gXwP5to2hGZC81hnVFVxKR/SRTabVerjOA6D6TjFYEpERDSBaIl2xPrXIh6vKtnWFF0wTQNuIKYcAqWUaE/YsIpToSV0JG0IKQDXHTzUVIphGAiHwzBNxpjxhv+PEhERTRQd/0DVHy9COG/BNA3lbqon7AFg5YYUfv1mn7dxSSAUNDB//nzlAMw9puMTgykREdFE0dcOKQFXD0HGGpW7JVuOUW67PWUDAGJBHRVBteCoSRcfmRodOPnP5fmJjMGUiIhoonAsAECy5gCkP3KZP48Y2DN67KxKfHxW6StJASCVSiEQCDCUEm9+IiIimjCcHAAJoftzkAkA7IGtpQEPdUkBcFmeAHDGlIiIaExzXRfOQEmnUvR8GkJISEOtHulQ2G7xilHv5Z+IGEyJiIjGsPXr1yOZTCq1bdz8LmpdF9IY/nqkRcVgGjC8BVPOmBLAYEpERDR2WWkYm59HOJtBOBwp2TyS2w5N0xCKxpH3aUjFPaYMpjQUDKZERERjlFx5A5rffhxSShiqd3rqGhBQL//0zo4c3u7MKbfvTDtwhYtMqh/d3eoF82tqapSfQeMXgykREdFYlewAAFixydBDlUpd3EAM6UmHKT/inr93I2urFcsvkkJicn0NZs9qVu4TiZSe8aXxj8GUiIhorLLzgJTo3O8MoHnhsL+9I+RgKF3WVqG8PB9wslg4vQENDQ3DPiYa3xhMiYiIxijp5CABSCMEP860552dV4T+ywHVyiftu7tdGNwzSkPAvzVERERjlTNwhMn0p/yT5RZmS03dW/knACz/REPCGVMiIqIykc/nsXnzZriuq9R+ajoBSAlpBJVnTPvzLjKW2p7R7kyhPmpQ9WDVAE3TGExpSBhMiYiIyoS9+R+wVz0LTdOhkuuEnYOmaTBDUahEzU29edzydCekLN12VyHTe8hk+ScaCgZTIiKicpDqRPjxKzA1m4VpGmp9BmYypal2on1Lvw0pCxWjwqZ6cFw6LabcFuBNTjR0DKZERETlIL0DkIV77DMNC5S75WvmQZpqdUmLxe8XNkdwziH1Sn2klOjt7UV3d7fymADuMaWhYTAlIiIqB06hiL0VqkP3wZf68oiBs0yeDjK5rgvDMNDQ0OCp1mgs5m2WlQhgMCUiIioPTh5SSkjTx3vsB2ZMDQ/BtLgs39TUxCL45DvuTCYiIioHA6WfXC3g2yPcgWDqYXvpYDDlYSYaCZwxJSIi8kkmk0E6nVZqG+zpRMh1IA3/ZkydwWDqfcaUe0ZpJDCYEhER+WHr6xBPfA9I9kJXCHXSzcBxXEDxINNQOEPYYyoHaksxmNJIYDAlIiLyw9rHEUhsAGxHvfyTaUCvne7bkJwh7DEFWDCfRg6DKRERkR+cPKQEepqOgmhbptRF6gFYVW3Kj9iatPH3zRkIxYr563sK+1iHsseUwZRGAoMpERGRH9xCCLRizXBr5/ryiIfe6MWarrynPlIKiFwG3d2OUvt8Po+amhoefqIRwWBKRETkB9culH/S/Ttl358vbBpd1BxBTURtu4Dm5HHU7FpMqqtWfk4wGOSMKY0IBlMiIiI/DJR/koZ/wdQaqJh/7KxKTK8JKfXp6+tDdUUELS0tvo2LaKgYTImIiBRZlgXXdZXaBqxs4US7j+Wf8k5hb2nIw6ZRKSUMQ/EwFtEIYzAlIiJSYOezWP+Pp2FZllL7lh1bEBQCMNRmMgEg5wj059SCL7AzmAYN76fsicoRgykREZEC/bfnY9qWVYVfq5ZbMnSEohXIKTTN2ALX/GkrsrbwPLaQh2DKGVMqZwymREREpQgXWs+6wq8DEUjFGUcn0oh81QyltjvSzmAoDZnqQXO/+jBiQW9L+TxhT+WKwZSIiKgU14aUAKTEux+9DUa4Ytgf4biFZfmGmIlrjmke9vffFZfyqVwxmBIREZUibACABACfTtnbQ7yVKZVKwXHUapIChQNcDKZUrhhMiYiISnHtwV9quj/fOl0xtINMuVwONTU1CIfDyn1isZinZxCNFAZTIiKiUoQDQEJqBuDTbOPOGVNv/TRNw+TJk1FdXT38gyIaYQymREQ0IfX392PLli1Kbc1MJ6ZaFoRm+rYM7gwEU9PDUr6UhT5cmqfxgsGUiIgmpGw2i87OTkSj0ZJtQ5k+CCGhB0L+BdOB8qUBj8FU0zQGUxo3GEyJiGhCklIiEAggHo+XbBtALwxDB0z1W5xcIbGhJw9r4LR9Ke8lCoX7vcyYAmAwpXGFwZSIiCamfBKBXDeMTOmC9ma2CwAgPRx8+tO6JB5+J+F5WAGPxfI1TWNdUho3GEyJiGji2bYKtY98GRX5LExT/RYkL8F0e7pwkr86bKAypBYcTUPDka3qJ+a5lE/jDYMpERFNPF2rAelAQoNQvsteQ7p5qfIjrIF77D8+uxLL2iqV+gghkEql0N9vKbV3HAem6d+BLKKRxmBKREQTj2NBSqCv4TBkDvuKL48o7i0Neaj/lMvloOs6YrGYUtgs7pMNBPwp+k800hhMiYhoXHAcR/kGJCOXhJQS0lA/zORVMZgGPdx7L6WEaZqYM2cOZ0FpQmIwJSKicWHdunVIJpNKbRvbN6HecSCVl/GBNV05rOnKK7ffkS6E5NAQDjMxlNJExWBKRETjQj6fhxBC6brNkCGhGzpCsThyCu8tpcSPX+xC3lEr/bSrWNDbiXmesKeJjMGUiIjGPLn9LTSs+TkgBALB0svz4b53oGma8oypLTAYSo9ojcFQrDVaHzXRWq2+XUBKyWBKExqDKRERjXny2TtQvfnvgKZB91Cg3g2WLq4PALa7s9bpp+bXKAdTr4pL+UQTFYMpERGNfVYKEkCy6XCgskmpiwhUINO8RKmtPXCQydDhWygFCsHUMNTrqhKNNwymREQ09jmFup/9U46CbDxg2N9+8IS9h4NMAGBZFixLrSYpAKTTadTW1np6BtF4wmBKRERjnnQLtyxJzZ9va8VgGvBQkxQAkskkKioqlJfn6+rqEI1GPY+PaLxgMCUiorKTy+XQ3t4O13WV2rekEoCUgGJd0m1JGz98fgeSeVG6MQBgIJgOYRm/sbERTU1q2wuIJjoGUyIiKjv5fB5dXV3KNxoJx4KhaTCCYahE2TVdefRm1ULvrmbUeivIr2kaT9kTecBgSkREZUeIwkxmPK52aj6gS2jQAUMtyDqiMAM6vymCTx5YrdRH14CqsPeDSTxlT6SOwZSIiMqOlB4K2UsJiIGrSHW1b2vFYBoL6qiN+vetkOWfiLzh+gIREZUdb8FUQBvYAyo9BlNzBL4LMpgSqeOMKRERlRc7i9gfLsaBO9YjYKp8m9oZYlVP5ReDqTECoZHBlEgdgykREZWX7vUw+jZCFzY0oXpqHrAqpkIqnsp3Bt424LEuaS6XG9z/qorBlEgdgykREfnOdV04jqPUVstlIIVAPlSPbYd/Xf0Z4RpAU1ubH5wx9VD+yXEcpNNpRCIR5T4VFRXKlQWIiMGUiIhGwJYtW9DZ2anUNtb7NqY5DkQkBDfa4Mt4HNd7XVIpJQKBAObMmYNwOKzcj+WiiNQxmBIRke9s20Yul0NVVVXJtkFDh6ZpCITVb0B6Y1sWf9uUglA8M7U1WbgpystFTsUT9oZhMGwS+YTBlIiIfCeEQDAYRDBYeg9o0AB0XVOuSQoAj63pR3uf+p30RbUR9W+DxWDKPaNE/mEwJSIi3wkhlAOdNlCT1Mu998W77I+bHcekCrV+saCOeY3qS/LFYMrZUiL/MJgSEZHvvARTr8XyC+9fCKb7N4Yxsy7kdXhKOGNK5D/+2EdERL5zXfV76TU5MGPqIZi6snjK3tu4vCgW/WcwJfIPZ0yJiMhfjoVI5ysIWllEUqVLLYX61gMApK6+x7RYl1T3EBoty0IqlVJun8vlUFtby2BK5CMGUyIi8sR1XWzatAm2bSuFtOq1/42mtb8BABgepjRVi+UDQ5sxzWQyqKysRGVlpXKfYDDIYErkIwZTIiLyxLZt9PT0wLZtpeLxFX0dAAA31gQ7Uqv0DGkEkWz5mPKY3IEZU69XjEajUbS2tnrqQ0T+YTAlIiJPinst4/G4UjANB3QYho5E23FITvu4L2MqzpiaHgvm84Q9UXnhf5FEROSJlHLwhLqKoZR/8koI70v5DKZE5Yf/RRIRkSdCCG/BdAin7L0aKGPqaSlfSgnDMHwaERENBYMpERF54n3GtHD9p5dT9l4IKSGLwdTjdzUeZCIqL9xjSkREnngNpsWC+V5mTN/ZkcOOtKP29sVUCu+HnxhMicqL8qfEN7/5TVx11VWIRqMAgN7eXtTU1Pg2MCIiGhlSSliW+j3zuVxuSHtMVW9y6kzZ+I/ndiiPp0hKgXSqH3nFA1Cq5a6IaOQoB9Prr78eF1544WAwbW1txauvvooZM2b4NjgiIvJfV1cX2tvbB0/bqxBCKLcdPPykuJSfyBVuiQqZGuY2qN9l3xy0UFVZgWBQrf5pTU3N4Pc0IioPysH0/R9YXj7AiIiofLmui1QqhdpatRqjVe89gVjPm9Db1TZ0mtlOAIDU1A4aFQ8y1cdMnHtovVIfKSV6enowdepUVFVVKfUhovLDPaZERBOclBKBQEBtplE4aFj3EDSoz5gWuWG17V/2QDINeKxJqmkal+aJxjjlYKppGpLJJMLh8OAHQCqVQn9//27t4vH4sA+SiIj842UFTBPOYCjt3v9s5QNNTnQSnOgkpba28F4sHwCDKdE44Gkpf7/99tvt9wcddNBuv9c0Da7rDu8IiYjIV14+t4ulnwAgNXUZoA1/1UHHHdotTpqmsWA+0RinHExXrlzp5ziIiGiUDOkWJ+i+hFIAcAZmTAMeat9zKZ9ofFAOpsuWLfNzHERENEpc1/Vck1S19BNQ2DO6pis3uHe0lE29hdJV3GNKNPF4PvyUSCTwxBNPYNOmTdA0DW1tbTj22GO5t5SIqEwkk0lPdUmz2ayv14v+YW0//rimv3TDXUgpYedz6O7uVmrvOA4qKioYTInGOE/B9IEHHsCFF164x4Gnqqoq3HXXXTjttNOGdXBEROTdpk2bkEgklO+BF0IgEokotfVakxQAejKFPnVRA1VhxfV54WL5zErMnt2q/Bxd1xEKhZTbE1H5UQ6mr7zyCs4++2ycccYZuOSSSzB37lxIKfHWW2/h9ttvx+c+9znMnTsXCxcu9HO8RET0IaSUcF0XsVjMl+Lxg8FUsSYpALgDe0aPnlGJ5TMqlfqkUikEAgE0NDR4HyQRjVnKwfSOO+7AKaecgvvuu2+31w8++GD813/9FzKZDL7//e/j3nvvHe4xEhGRIiGEt3vsPdKGcO+9M1Dy1Ospe56wJ5p4lD9ZnnnmGfzoRz/a69e/9KUv4fzzzx+WQRER0dBIKT2FOj2fQM2aX0G3U0rtDSs90FE9mBZnTA2PdUkZTIkmHuVPlo6Ojt3qmL7ffvvthy1btgzLoIiIaGiKwVR1xjS2/UVUdDzt+TlOSO0WJwBwZLEuqfr7SymV98gS0fihHEwzmQzC4fBevx4KhZDL5YZlUERENDRCeLsqVLOzAIBczVykJh+p2ElDrn6+8jPcIS7l84Q90cTj6VT+H//4R1RVVX3g1/r6+oZjPERE9D5eSj/l83kIIZRDne7mC8+onIb01KOGNL5SnCFcMco9pkQTk6dgetZZZ33o1/nTLRHR8Eomk1i/fr2na0Mty1I+ka8NBFNp+FdmyRESrnCRSvajO5BR6mNZFoMp0QSkHEy9Lg8REdG+cxwHmUwGFRUVyn1CoRBMU+3jvThjKnwNpoAUErXVVZgxY5JyPz/KXRFRefN88xMREY0cOXBwSLVwfHzDI4hvegyQapMJ+uCMaVB5TPe/0o03t6ufKcgO1Iuqra5CU1OTcj8imniU10lefvllHH300Xvc+gQUrik9+uij8dprrw3r4IiIJrpiMFVVseWvMOwkDCet9I8mHUhosOLTld7fdiVe2pxBxhbK/0gJBAygKc5bmYjowynPmN5yyy342Mc+hng8vsfXqqqq8PGPfxw33XQTHnjggWEdIBHRROb1dLruFGYydyy8CHbFZKU+biAGEapWamuLnUH5qmWTlA80OZkEqqPq15gS0cSkHExfeOEFXHXVVXv9+kknnYT//M//HJZBERFRgefyT4On7FvgxIZ/2dzZJZhOiQeUQ3OPrfOALBGVpLyUv2XLFlRW7v2O44qKCmzdunVYBkVERAXFgvmKjXeesjf3Xnd6XzjuzmL5XoImyz8RkQrlGdOGhgasXr0abW1tH/j1d955B/X19cM2MCKi8SibzSKbzao1lhLuqv9FU+dGxLpiJZtrUkBDYYZVGD4F0+L1oprmfTaXM6ZEVIJyMD322GNx/fXX4/jjj9/ja1JKXH/99Tj22GOHdXBEROPN1q1b0dHRoTR7GE5vxn6v/wRVmgbDVL+eU+ghT6fsvXAlICEhXQeJREK5XygU4owpEZWkHEy//vWvY/HixViyZAkuu+wyzJkzB0BhpvSWW27BmjVrcN999/k1TiKiccF1XQSDwb3eorersOyAGTDhBuLon3y48jNydQcAmj8h0C4kUwQCBmbNmvWhV1W/XyxWetaXiCY25WA6c+ZM/OlPf8LnP/95fPrTnx5ckpFSYv/998cTTzyBWbNm+TZQIqLxwMt1oZqwAQBOtBG9cz/j57CUFZbyJUxdRzweRzDoz8wsEU1MngrsH3LIIVi1ahVeffVVrF27FlJK7Lfffli0aJFPwyMiGl+87MvceV2of2WWknkXv3s7gbSlNq60JSDh/fATEZGKId38tGjRIoZRIqIh8DRj6hZmTKXu36zkPzqyeK497blfddhgMCWiYccrSYmIRpDrusqBThcDM6a6fzOmxetCZ9QGcdhUtT2gtm1hTp16DVMiIlUMpkREI6XvPdSv/jngWggESofNYGozAEAY/l3laQ/UJZ0SD+LI6RVKfTKZDDRN4yl7Ihp2DKZEREPkui46Ojrguq5S+6pX70L15icBALqhHurcUOkT/ENlu3Kg/JONdFptST+TyaCyspIzpkQ07BhMiYiGyLIsbNu2DblcTmkGNNS7HQEpkW86BFb1TKVnSCOI9OQj9nWoe2ULCeEKBA31GdCKigqWfiIiXygF09dff135DRcsWDDkwRARjSXF60Jra2thmqU/TsMBDaZpoK95CdIe6pL6yXIL5Z/iFVEsXLhwtIdDRBOcUjBdtGgRNE3b633Nxa9pmqa8pEVENNYVg6nyYaaB8k9+7hl9eUsGT29K4YM/rfe0PWVDSiAcMLk0T0SjTimYbty40e9xEBGNOcVgqmpnXVL/gulja/qxLWl77tdQ4d/JfyIiVUrBtLW11e9xEBGNOeaLP8LsNX9FIKBWZzSQ2Q7A32BquYXyT/88rwr1MbVjBE42iYOmVvo2JiIiVUM+/PTWW2+hvb0dlmXt9vrJJ5+8z4MiIip7+STMd36HUC4HwzCgugouNQN2tMG3YVkD5Z8OmBTB5LjaLGh3dxamYfg2JiIiVZ6D6YYNG/Av//IveOONN3bbd1rcm8Q9pkQ0lil/hmUTAAChmeg69Erl93eijRCh6iGMTE2xLmnA8LZflPtLiagceA6mF198Mdra2vDnP/8ZbW1tePHFF9Hd3Y3LLrsMN998sx9jJCIaEV1dXdi8ebPSvtFQ6j1Mz+fhGjHk6/YfgdGVJqWE5Rb2vQo7j1zOUepn2zaDKRGVBc/B9LnnnsOTTz6J+vp66LoOXddx5JFH4oYbbsBXvvIV/OMf//BjnERE/sr1I/CP+xDv2IhwOFyyuWn1A5oGPVw+ezNdCUhZmPXVIaA6+RuPx5XKXRER+c3zJ5HruqisLHwQ19fXo6OjA3PmzEFraytWr1497AMkIhoR6/6EyOpfo9GyYJqK+y11DTJa59uQbFfi3pe70ZVWm/ncdZ53Vlsr6mtrlJ/FYEpE5cDzJ9GBBx6I1157DW1tbViyZAluvPFGBINB3H333ZgxY4YfYyQi8p+VAgBkYtMgWpYodZHQkWk6zLchberN441tWY+9JKrDOqLhEMMmEY05nj+1vv71rw/ep/ytb30Ln/jEJ/DRj34UdXV1+NWvfjXsAyQiGhGuBUAiWzkduZmnjPZoAOw8Yd8YM3HaArXZTyEE4lpW+XpRIqJy4jmYHnfccYO/njVrFt555x309PSgpqaGm+eJaOxyCzcgSb18Cs0XT9jHQjrmNJTe9woAjuMglbL4eUxEY5LnH6kTiQR6enp2e622tha9vb3o7+8ftoEREY2ogRlTqZfP8rctCsE06LH0k67rDKZENCZ5/gT+9Kc/jZNOOgnnn3/+bq8/+OCD+N3vfodHH310SAPJ5/PI5/ODv2fIJaJ9lclk4DhqB4eCmSSkkGU5Y+olmBZLXXEpn4jGIs/B9IUXXsCtt966x+vLly/H1VdfPeSB3HDDDbjuuuuG3J+IaFdSSqxfvx6pVEqp/dRtHahxXWimf9eFelWsSeraFvr6+tT6WBYikQhnTIloTPIcTPP5/AfOQNi2jWzW6+nRnb761a/i0ksvHfx9f38/Wlpahvx+RDSxSSGA1HZUSIFIJFKyfcxwYRgGQtFKWCVbD80b27L471V9cETpAv4AkHMEpJQIBww0NKhfY2qaJgKB8pn5JSJS5TmYHnbYYbj77rtxxx137Pb6XXfdhcWLFw95IKFQCKFQ+cxUENEY9+S3MXPV7wtF8HXF2UMNvu4xfWlLBt0Zta0FO0lMrQ6xHB8RTQieP4G/853v4Nhjj8Vrr72GY445BgDw5z//GS+99BIef/zxYR8gEdGQbH+zUHBeNyEUw6YIxZGrnefbkCxHAABO2C+Ohc2lZ3EBwMplMK2uwrcxERGVE8/B9IgjjsBzzz2Hm266CQ8++CAikQgWLFiAe+65B7Nnz/ZjjERE3gkbAPDe4q9Bqy+Pz6ZiXdJJlSamVgWV+vRrOR5kIqIJY0hrVosWLcLPf/7z4R4LEdHwcQeWzMuo/JM1eMpePWhKKRlMiWjCUPrE7u/vRzweH/z1hym2IyIaVW5hxlQzyucQkDXE8k+GYfg1JCKisqIUTGtqarB161Y0Njaiurr6A8uQSCmhaRpc1x32QRIRWZaF7u7uwTqdpdTmClcnSyMIvwon/XFtP9Z05ZTb70gXZnGHUjCfiGgiUAqmTz75JGprawEAK1eu9HVAREQfJJ1OY+PGjWr1OaVE3CoERt0MQi3KepN3BB5+O+G5n3BdyGwC3d1qs6Cu63LGlIgmDKVgumzZMgCFO5ifeuopfOELX8DUqVN9HRgR0a6KM6XFH5I/lHBgmoUwpxmmL8G0eCsTAJx1cC1U69kH7DRmTmlEXV2d8rOi0ajX4RERjUmeTgWYpombbroJZ555pl/jISL6QMXtQio0sbNWqF9XjDoDuVTXgEOnxpT7dXfnUFtbi/r6el/GRUQ0lnk+rvqxj30MTz31FKZPn+7DcIiIPpi0swjke2FkS89/6k5mZz/Nn1P57sDtTaZq8f5d8LpQIqIP5vkT+4QTTsBVV12FN954A4sXL0YstvtMwcknnzxsgyMiAgCkdiD+27MwN50YXKJXIaEBmj8Hh4pL+V6DqaZpDKZERHvhOZief/75AIBbb711j6/xVD4R+aJ3IzSncJjJy9J8etJhUN786dHOGVNv/ViXlIho7zwHUyGEH+MgIto71wIgkamYjt5l14/2aADs3GNqcCmfiGjYlM+VKEQ0oSSTSWSzWaW2ge5ORFwXIlg+ZZOcgaV8QwNs21bqU6wswGBKRPTBlILpD37wA5x33nkIh8P4wQ9+8KFtv/KVrwzLwIhofNu2bRs6OjoQDJa+M756x3tocVzoZsjXMWVsActRWxXqy7mQUkK4tnLABgqln1iXlIjogykF09tuuw1nnHEGwuEwbrvttr220zSNwZSIlAghEIlEUFVVVbJtRSYM0zRghv2r5/nOjhx+9PwOCE9FTyVCpoE5c+YoBWyg8Dmp2paIaKJRCqYbN278wF8TEQ2VEMJDXdLCUrnU/dt99G6vBSELZ6VUbwyVUsP8xiAikQgCAX/qpRIRTST79CnP/VJENFSegqlbKJgvdf9mGq2BPaNHTa/AJ+fXKPXJ5XJwXZefgUREw2RIwfSee+7BbbfdhrVr1wIAZs+ejRUrVuCLX/zisA6OiMYp4SLYsxq6lUXYKr08H0h3FH7h44yp5Rb2lgZVp0ux8zYqBlMiouHh+VP+m9/8Jm699VZcdNFFWLp0KQDgueeewyWXXIL29nZ861vfGvZBEtE48+LdmPLyTyGlhGGo1/T063pRYOeMadD0HkxZl5SIaHh4DqZ33nknfvKTn+D0008ffO3kk0/GggULcNFFFzGYElFJsu89AIAbqoYbjqv1MUJITT7CtzENBlMvQXmgWD5nTImIhofnYGrbNg455JA9Xl+8eDEcxxmWQRHR2OI4Dt59913Ytq0U0hq6tiMkJXbMOBVO29G+jOnhdxL4y4akcnvLlZBSwrWy6O9XO5qfyWRQU6O2H5WIiErzHEw/97nP4c4779zjStK7774bZ5xxxrANjIjGDsuy0N3dDdd1YZqlP1Zq8hmEAARCYfj14+yL76WRdzzVfoKUAtNrI6ioiCm1j8ViqKysHMrwiIjoAwz58NPjjz+Oj3zkIwCAF154Ae3t7TjzzDNx6aWXDrZ7f3glovFJysJsYzweVwqmIVMr7C31sWB+biCUXri0AfVRtY+6XCqB2dMbMW3aNN/GRUREe+c5mK5atQoHH3wwAGD9+vUAgPr6etTX12PVqlWD7bjnimjiKAZT73VJ/TnMJKVEbuAGp6bKAKrDajctdef42UVENJo8B9OVK1f6MQ4iGsOEKIRA9bqk3oPpu30WejJqC/+OkBgos4ywh1P2AHjCnohoFPlXFJCIJgxpZxFMbUbQSCqFU83NFfoZasF0a9LGTX/d7nlcugaEPNQlBRhMiYhGk+dgmsvlcMcdd2DlypXo7OwcnCkpeuWVV4ZtcEQ0BrgOog+fh/16tsA01ZbMi6Sm9hHUPTBTGjI1tFSp3/60oCnieWmeS/lERKPHczA955xz8Pjjj+Pf/u3fcNhhh/FDnGgcklLCsiy1xpku6OkdAAA3WKX8DLtiCuyKKWptB2qMTq0KYsURjcrPGAp+phERjR7PwfSRRx7Bo48+iiOO8K/QNRGNrs7OTmzevBlSli63FMh2YpZlw9WD2Hb0Hb6MpxhMA7q30Njb26v071AkhGAwJSIaRZ6D6ZQpU1i3j2icc10XqVQKtbW1JdsGrcIsoxFSq/05FMVbmQIe77GXUqKhoQEVFRVKfTRNQzyudhMVERENP8/B9JZbbsGVV16Ju+66C62trX6MiYhGmZQSgUAAwWDp/ZxBQ0LXNQgfa5IONZhqmoZJkybxh2kiojHCczA95JBDkMvlMGPGDESjUQQCu5+q7enpGbbBEdHoEK4DDQIQbsm2ulM4YS909UNJaUvgje1ZOK7aMvv6njwAIDiEYMqleSKiscNzMD399NOxZcsWfPe738WkSZP4oU803mx9HQ2/vxg12aSnU/bSUJ8xffidPjy9Ke15aF5qkhaDKcs/ERGNHZ6D6bPPPovnnnsOCxcu9GM8RDTatrwMzcl67parO0C5bV+2MBM7tSqA2ojax1DI1HDUdLW9ogBnTImIxiLPwXTu3LnIZr1/0yKiMcK1ICXQ0/RR5Bd8VqmLhA4ZiCo/orhn9NhZlThkij+HphhMiYjGHs/B9Hvf+x4uu+wyXH/99Zg/f/4ee0x5opWo/Agh9rgMY280KwMpJUQgBhFQn6H0ohhMg4b6MrvrushkMurPsCyEQiEu5RMRjSGeg+nxxx8PADjmmGN2e704O+G6pQ9LENHIkVJizZo1yisdk7a0o0oISEP9MJNX9mAwVZ/NzGaz0DQNkUhEqX0kEkEoFIJp8uZlIqKxwvMn9sqVK/0YBxH5RAiBbDaLfD6vFOoMYUPTNASjcSje/eSZPcTyT5FIBPvvv79PoyIiotHmOZguW7bMj3EQkU9k9wbUvPsYdF3fY+vNB4lkO6DrGmCG1d5fSvzitV6s684rj6kr4wDwdpOTlJLL8kRE49yQ1rj+9re/4cc//jE2bNiAhx56CFOmTMHPfvYztLW14cgjjxzuMRLRPtD/8h00bHlroHSSehAUAbVDSWlL4Ll276WfAoaGuqh6OSoADKZEROOc52D6P//zP/jc5z6HM844A6+88gry+cIsSSKRwHe/+108+uijwz5IItoHmV4AQLrhIGhhtcOJbrAK2Qa1knC2KCzLGzrwlcMblYdVHzURC3qok8oZUyKicc9zMP3Od76Du+66C2eeeSZ++ctfDr5+xBFH4Dvf+c6wDo6IhoFT+OGxe9YnoVW3DPvbD+4X1TXMrPXvWlIpJQzD2wwrERGNLZ6nH1avXo2jjjpqj9erqqrQ19c3HGMiouHkFo4webmZyYvijKmXg0xDUaz8QURE45fnGdOmpiasW7cO06dP3+31p59+GjNmzBiucRHRXriui1QqpdZYCkTsPCAlYJQ++DQUzkB5VNPD/lUASCaTsCz1c/+WZXHGlIhonPMcTM8991xcfPHFuPfee6FpGjo6OvDcc8/h8ssvxze+8Q0/xkhEu+jv78e6detg23bJtpprYb5lFWYbTX9mTJ1dlvK9sCwLdXV1qKqqUu5TUeFPwX8iIioPnoPpVVddBSEEjjnmGGQyGRx11FEIhUK4/PLLcdFFF/kxRiLahRACjuOgrq6uZFvdSsE0C7OMfgXT4lK+6XEpX9M01NXVob6+3o9hERHRGOQ5mGqahquvvhpXXHEF1q1bh1Qqhf33358zGUQjwc4h8uxNmLF1HcIhhTqjslAvVGoGoKstg3elHTz0Ri8yjtoVplm70M7rUj73jBIR0fsN+a6+YDDIG1iIRtrW1xBo/xti+TzMnPp+SyfSoNz25Y4M3uzMeR5avceapAAYTImIaDe8RJpoLHEKgTEXmYTMAWcod7OqZ6q3dQpL8wdOCmPpNLUi+7qmYXa9960CDKZERLQrBlOisWSg9JMTrEZ20mJfHlHcM9pUGcDC5qgvzyhiwXwiItoVgynRWOLakFJCaP6VTbKHeMpeCLU9qcDO/aWcMSUiol0xmBKNMsuylOt5GukENCEgdX9qkgJDK5hv2zb6+vpgmuofKeFwmDOmRES0GwZTolH23nvvYceOHUpt6zo2otlx/A2mAzOmXk7ZCyEQiUTQ1taGYDCo1EfTNESj/m4VICKisYXBlGiU2bYNTdNQWVlZsm2024RhGAhFK5BWfH8pJTK2VB5PbqBMlJcZ0+LSfDwe9zRrSkREtCt+ByEaTY6FUOdrqMqmEbUiJZuH01ugaQAMtVlJAPivf/Tgpc0Zz0PzsseUe0aJiGg4MJgSjabnf4imf/w/CCFhGOr7LaWuHkzf2eG9JmksqGNGrfozGEyJiGg4MJgSjab+rZASsMP1cCLVSl2kGUZqyhHKj7AG9ox+42NNqI+q/SevaYXapKqklNB1nYeZiIhonzCYEo0iOVCXdEfbKRBty4f//aUcDKYRU4fhsQSUl+dwtpSIiPYVgynRaHItQErA8OeUvS0Kbw8AQVM9OGazWWQy6vtSLctCfX291+ERERHthsGUaBgJIdDV1aVcbD6e6gfgbc/oloSF3pyr1Dbv7DyNH/Rwyj6bzaK2thZVVVXKfcLhsHJbIiKiD8JgSjSMcrkc3n33XeTzeRhG6duZZif7EJYSRjAMR+H9t/RbuOGp7Z7HFTQ0z3tGY7EYmpubPT+LiIhoqBhMiYaREAJCCNTW1ioF05CpwTQNGAG1YNqTKcyUBk0NTRXqy/8HTy5diur9eJCJiIhGGoMp0TCSUno6CKSJQhyVintMnYHrQqdVBbHiiMahDVIRgykREY00BlOi4SIlAs/djlnrX1C+ltPI9xW6Kl4xWgymHkqeDhlP2RMR0UhjMCUaLpkemOv+iHAuB8MyoJrrpB6EE6pWausOnKkyRiA0MpgSEdFIYzAlGi5uHgAgNBNdh16p3M2JNUEGYmqPkMUZU2+hsbjFwEt7LuUTEdFIYzAlGi6uDaAwA5qv29+XRxSX8k0PmVFKia6uLpim+n/uwWCQwZSIiEYcgynRh8jlckilUkpt9d5ORB0HQvd+Al5VcSnf9DBjKqWEaZqYNm0aKisrlftFo1GvwyMiItonDKZEH6KnpwcbNmxQKv0USb2LWY4DhPy5xQkAXOF9Kb9YJSAej6OiosKvoREREe0zBlOiD1Hca1lbW1uybUjrgmkakCH/Zkyd4h5TD1tMi8GUS/NERFTuGEyJPoQQwkNN0oE9plrp2dWi7oyDF95LDy7Rl7Kuu3DAaigzpjxlT0RE5Y7BlOhDuK6rHujEwP31ijVJAeCRdxJ4aXPG87giAfXZTwZTIiIaKxhMiT6E1r8FkdRmBPXekm2DqS0AAKmr/2eVsgpTpXMbwphUodYvZGpYNl19ryiX8omIaKxgMCXam7V/wqS/fB11jgPTVF+e9xJMHbewZ3TptBgWT1E7BW9ZFpLJPnSn1Z5hWRbi8ThnTImIqOwxmBLtTe9GAIAwwnDCcaUuUjOQmrpM+RGDh5k8TGbm83lUVlaivr5euY9hGJ7qmBIREY0GfqeiCSWZTMJxHKW2oWQfNCnQ23wUcgvP8mU8xRnTwBDqkjY3N/syJiIiotHCYEoThuu62LhxI5LJpNKy9pRtW1DruNADId/GZA+cxvd6xSj3ixIR0XjEYEoThpQSjuOgoqIC4XC4ZPvY1iACARNmyL8bkIoF873OmDKYEhHReMTvbjRhSCkHT6grGahL6qX8k1fOQDA1PVTMl1Iq3URFREQ01nDGlCYMOXDQSNVgwXxD/T+TvqyDrUm1PawAkHcGgukQbnIiIiIabxhMacIQQngKdYPBVHHG1HIEvrNyG3KOtwAMFK4YFULt+ichBGdMiYhoXGIwpQlDW/tHTHv9IRimqbRHM9i/CYB6ME3bYjCUTqlSX/6vCzgIOikkEmqBWdd17jElIqJxicGUJgzzlXsR69sC3TDgZSXcCdeptRuY8AyZGr66rEn5/bu7u9HU1IS6OrXnAEAkElFuS0RENFYwmNKY5bouMpmM8t7RUDYFCaB7v09DRmrUnhGqRr5mjlrbgYNMxhD2f8bjccTjakX8iYiIxisGUxqzEokE1q9fD9d1ldofkE0BUiLTdBhkVP3WJFWDJ+yHsMrOw0xEREQMpjSGCSFgWRZqa2tLN5aiEBh1E5oZhvfjSaUNXOLkuVg+wIL5REREAIMpjWHFE/YqoU5zrMF9pVL356/9zhlTbzVJAc6YEhERASywT2OYl7qkmthZW1QaQT+Gs3OP6RBucWIwJSIi4owpjWFSuNCEBc21SrbV7XShj2YCmj8/jzlCQkLCtS10d3cr9RFCwFQsX0VERDTeMZjS2JTcjuqHz8b8ZBdMU73YvF/L+ADgisIMaMA0MHXqVJim2rM0TUM4HPZtXERERGMFgymNTTvehpbv99wtWz9fuW3OFni2PY2co3Yj07ZUYbuAqeuYPHmycjAlIiKiAn7npLJh2zZs21Zqq2eSgBBIxWehf+lXlZ8hjZBy22fa0/jfN/uU2xceAIRNtQNZREREtDsGUyob7777Lnp6epTaVm9fjymOA2mEIE1/lsETuUJ91ClVAbTVqB2YcmwbS1uiPMxEREQ0BAymVDYsy4IQApWVlSXbRoIGdMNAMFLh33jcwhL+wqYITpxTpdQnnU7DNE0GUyIioiFgMKWyIaWEaZpKezMNiEJdUh8PM+WdQvmngOG9/BMRERF5x2BKZcN1XeWZxmJdUqkHfBuP7UpIKSFtC6lUSqlPOp3mCXsiIqIhYjClslG8yUmFJgqHpLyUf3rxvTRe2pJRbt/eZ8F1XUSCBgIBtQBcXV2ttBWBiIiI9sRgSmVBCoFAqgMB4SAgSs84GrnCISkvM6a/eTuB/oEDTV60NddhwYLZnvsRERGRNwymVBbk8z/C9BfvBzQNupcrPT0E02I90n89oBoVQbV9oDKfwtzGqPIziIiIaOgYTMk3lmUp32evd64BAAgjAhlQqzUqzAgykxYrtZVSwnYLYzlkahTxkNptUd3dOR5mIiIiGiEMpuSLbDaLtWvXwrJK32MPAK3dOxCREp3zzoLdcsSwj8cWQDEjBz2csgfAYEpERDRCGEzJF67rIpPJIBKJKAU7U5cwDB1mKAK1u5+8sXa5VtRLMNU0jTVJiYiIRgiDKfmiuIQfDAaVgmmhLqkGqan9lXSFxMoNycHbmUop1iQ19cJBK9UjUF5KWBEREdG+YTAlX0gpPZV/wkBdUtWC+Wu78/jNWwmvo0IQAslkUrlHOByGYajtRyUiIqJ9w2BKvhBCeKtLKgtzmKp1SdNWYWm+JmLg0Klqp+aFkJgeczB9+nTEYjG1cWkaIpGIUlsiIiLaNwym5AtPs6XY9SYntdlJVxSW5idVBHDyvGqlPo7jIJ1Oo6KiQjmYEhER0cjhcWPyRXEpX9VgMFXdYzrw1oaHv8HFsMxT9kREROWJM6akREoJ11W/Ncm2PZ6tl972mDqieJjJQzH+gaDMw0xERETlicGUlHR2dqKjo0NtFlRKTFr3S8zqW49QSK1Yvm4X7rBX3WNaXMo3PAZTXdcZTImIiMoUgykpsSwLiUQCNTU1Jdua2S7UbH0KmgZoeS/XiwbhBuNKbQeX8oeQMbmUT0REVJ4YTEmJlBKhUAjhcLhk24AtoesahBnFjvnnKT/DrpgCaaqdgHeEhJQCVi6L7u5upT6WZaGqqoozpkRERGWKwZSUCCHUSz+5+UKfQAzZxoN9GY8rJISUiEbCmDFjhnI/0zRhmvxrT0REVI74HZqUeAmmejGY6kHfxuMKABIImgaampp8ew4RERGNHAZTUiJcF5Ci8E8JmpMDAEiz9LJ/0Zvbs3h8bRKuYompnmyhQkDAS70oIiIiKmsMplTaey9h8uOXozmfhmGqX88pDbUT+QCwckMK63vynofWUOHfrCwRERGNLAZTKm3Ly9DdHFyPh4aytfsrt805hZnY42bHMb1GLWzmMyl8ZKbaKX4iIiIqfwymE5CUEl1dXcoF86OJXphComfy0cjtf5raMzQdMqB2hz0AWAP1n2bXhzC3QW0LQE9PFgEeZCIiIho3+F19AnJdF5s3b0YikUAwWHp2ckrXdtQJAS1UCRGs8GVMeUdASAEnn0EyqXZrVDabhWGoby0gIiKi8sZgOgFJKSGEQG1trdLNTNGtQZimASOofphpU28eHf3q15KmLQHhClREQohG1Z4TiUSU6qoSERHR2MBgOgEpXSu6C00OLPnrarOTaUvgtmc6CyWdPGqb2oz9pjV770hERERjHoPpBCSlhJRSvWC+cAr9NLVgmrJcuALQNeCASeozmnWmjYZKzoASERFNVAymE5AQYmjBVFf762IPHGSqCOr498MalMfV3d3N60KJiIgmMFYnn4C8zpiiuJSvOGNaDKYBQz1kFrcXMJgSERFNXAymE5DnpfyBYKo8YyqKwVT9r5eUErquM5gSERFNYFzKHwcsy0IikVBun06nPR2AGupSvqFJ2LbayXwhCieldJ0/KxEREU1UDKbjQDKZxNq1a5XbN27+I2Z3vYRgMACg9AylmeuGkEDO1ZC2ShflT1mFPaxwbWSzWeVxRSIR1iUlIiKawBhMxwEhBDRNQ21trVL7yf94CrqbBRy198/YEl1pB9/9u4YteodSHyklwoEA5s6dC1PxdiZN01iXlIiIaAJjMB0HintGFRtDc3MAgM6DVkCYsZJdVm5I4vHNJrr0OuUxaRqwf0MQsViMy/NERESkhMF0HPBUMF840FBon6udB2lGSnbp6OhFl57Cx2dV4qR5VUqPSafTMA2Dh5mIiIhIGYPpOOAlmOpufmc/o/R1pAAGb3AKGBp01ZP8AE/ZExERkScMpuOAdPIwnDR0u3TQNPK9hT56ANDUltgdUTxl7y1kcgmfiIiIvGAwHev63kPtI2ejMtMP01Q/0S70oHLbwWDqIWcW65ISERERqWIwLUNSysG6niV1vg3Y6iWZirKTFiu3dSUghItcJo2eHrWj/LZto7Ky0vO4iIiIaOJiMC1DGzduVC6YH9+2Fs2Og2T1/uhfepX6Q3QPs6tCQgiJqnglZs5sUe7H0k9ERETkBYNpGcpms0in06ioqCjZVpcuoGkIhqOewqYXxaX8ylgUjY2NvjyDiIiIiMG0DEkpEYlEEImULuUUMnXougaYaifsh8IZPJXPPaNERETkHwbTMqS8vxSAJgp30UtNfbY0lXexqc9Sbt+fL1xDahgs/URERET+YTAtM8VbnJTrf4rCYSSpq/9f+YPndqCj3/Y8tpCHU/9EREREXjGYlhmvwVSTA6fkPQTTrnShz5R4AIau9pwwgAObS+95JSIiIhoqBtMy4zmYFpfyPQRTW0g4joNPzw4hHlabBdX1IKKhgPIziIiIiLxiMC0zMtWJ6s1/hq5JBAKli+CH+tYX+ulqodEVEmIg+LZOm4rGmrjy2OJx9bZEREREXjGY+sx1XXR3dyvfZx99+v+isf0ZaJpWOG2vSBhqNUNtIYGBoTQ11qMmziL4REREVB4YTH2WzWaxadMm2LatdEXnzO7NiEoJq24e3Eid0jOkGUF66lFKbR13Z0AO8jATERERlREGU58Vrxetra1VCqbhgAbTNNA98yTk6hcM+3h23nuvwTQYTImIiKh8MJj6TAjh7TCTW6gvKvXS+0uLujMOkgO1Rkvpy7kAJEwd6iWpiIiIiEYAg6nPhnzK3lALpu19Fm7863ZvYwIQMDQGUyIiIiorDKY+Uz30VDQ4Y2qoXTHakbQhpYQOgYqgWtAUQuDwlgqlrQVEREREI4XB1GfGu39DU/vziPXElNrrThaAevmnvCMghMCC5ggu/OhU5XEFAgEGUyIiIiorDKZ+Snch+uz/RWMuD9PDCXgJDSIQVWqbdwozshXhAGbMmDGkYRIRERGVAwZTj6Qs3JqkJNUNKSSEHkD/tGOU339roBWd6SAAq2T7zoHrRUMBnrAnIiKisY3B1KNt27aho6NDqW04uQnTbRtuoAq98z6n1OeXr/fg6XVp4G0vB5okokEGUyIiIhrbGEw9sm0bqVQKNTU1JdsGBkoyGSG1ZXkA2NJfOJUfC+oIGmqHmQwpcchU3uBEREREYxuDqUdCCIRCIYRCpU/Nh8yBa0UNtYNMAGAN7Bn9/MF1mNeods1od3c3Wuoiys8gIiIiKkcMph4JIZTbDtYkVTxhDwB5V8IVLnLpfvT0ZNT65PMweIsTERERjXEMph4JITzc4lQMpup/zJYjIIVEfU0Vpk4uvV2gqLKSS/lEREQ0tjGYeiQdC7qwoDm5km01J4ucLbElpeHpDUml98/YhaX8yU2NaGlp2qexEhEREY0lDKZebPwbpj5xFSY7FgyjdHF6VwBb0g7WZB38d7rP06MqQvy/hoiIiCYWph8vOl4BhGINUwCukBDQsNqcg4OnKJ7MlxINQRtVEfV9qURERETjwYQPprlcDq7rKrUNZNOQUmLHtBNhzTu1ZPuNPXnc/mwXqmIRXLe4TukZQggkEgleF0pEREQTzoQOpq7rYs2aNchms0rtp2zdjGohgEAU0ihdLioHCUcLKNcjBQo3PwFQPmBFRERENF5M+GBq2zaCwSCCwWDJ9mETMAwDgUgF8grvb7uyEDRdC93d3Uph03VdhEIhBlMiIiKacCZ0MJWyEBxN04Rplv6j0KUDy5XYnAK27Sh9Kn9TrwUpJcLBAKZNm6a8PK/rOsJhteL6REREROMFg+nA0rmKrX0ZBFIOfvNOBn9ft0O5X8jQMGXKFM6CEhEREX2ICR1MtXceRusbv4VhmkqzmYnUu7AAmMEQplSqnZqXrosjpkUZSomIiIhKmNDB1Pz7TxDr74FuGFDJjRnhwgKw37Qp+OwCteL3qVQKoVDpg1JEREREE92EDqZw8pAAuuacAS1U+krPP63rxwvdESyItSo/QkrJ0k9ERERECiZuMBUCkAIAkG5eCi0cL9nlzS1dWGdksFAWTvOrsCwLFRUV+zRUIiIioolgXAVTx3GQSqUUG+cQdQducdLV/hhsV8J1XQjHVq59GggElEpREREREU104yqY9vX1Yf369Uo3ORlOGgdYNiAlYKgdZLKFhJRAfW015s+fpTyuQIDXixIRERGVMq6CqZSFGc26utLXf+p5E7bUkbUF/vvNJFROP21LFWZYY+EQ64wSERERDbNxF0xVyzJpwkZ3xkFemnhqU9rTc6oinAElIiIiGm7jLpgqF8x3bQgJuNBx9IwKBA21k/Omk8Hsxug+jJKIiIiIPsj4CqZWBoF8H4xc6VlTLdsDAHBg4oQ5VYgG1IJpT48L0zD2aZxEREREtKfxE0y716P2kS+gIpuGaZYOjkICmwHYMGAlepHT1bYAuK7LW5yIiIiIfDB+gumO1YBbqC0qNZVgKuFCx8v6AnysrhaRSETpMZqmIRrlUj4RERHRcBs/wdS1AACJmgVIHX5lyeY9GQffeKIDOgTObm1VDqZERERE5I/xc1emawOQkIo1SR1ROCQVMDReGUpERERUBsp+xnTjxo2orCx9j73+XjuQyOP5lMBv/tRRsr0jAEDC1DXuGSUiIiIqA2UfTNvb29HY2FiyXXZrN5odiTRMdGdK3/wEABJAU4XJYEpERERUBso+mFZWVqK6urpkO9vQoWkaamJhXHJovdJ75/M5TIkHuZRPREREVAbKPpiqkq5duPlJ11Fn5tU6mRqiYQZTIiIionJQ9sG0ee3PUbej9B5TLbkaEkAgFMEBBxyg/P6myaV8IiIionJQ9sFU2/IS3N7SJ+1DtkAOgAjXoKKiwv+BEREREdGwKvtg+oB7LKKuQo1RHUgZQbRNOtL/QRERERHRsCv7YJpsORpaVa1SW83OYumsST6PiIiIiIj8UPbB9LhpBpqa1G5l6u3NoToa9HlEREREROSHsg+msVgMUkqlttXV1QgE1G5+IiIiIqLyUvbB9MADD0Q8HldubxiGj6MhIiIiIr+UfTA1DINhk4iIiGgCKJtgms/nkc/vLIyfSCQAAP39/aM1JCIiIiL6EMWcprrtspSyCaY33HADrrvuuj1eb2lpGYXREBEREZGq7u5uVFVV7fP7aHK4Iu4+ev+MaV9fH1pbW9He3j4s/6JERKOtv78fLS0teO+99zztnSciKleJRALTpk1Db28vqqur9/n9ymbGNBQKIRQK7fF6VVUVP8CJaFyJx+P8XCOicUXX9eF5n2F5FyIiIiKifcRgSkRERERloWyDaSgUwjXXXPOBy/tERGMRP9eIaLwZ7s+1sjn8REREREQTW9nOmBIRERHRxMJgSkRERERlgcGUiIiIiMoCgykRERERlYWyDaY//OEPMX36dITDYSxZsgQvvvjiaA+JiEjJDTfcgEMPPRSVlZVobGzEKaecgtWrV+/WZvny5dA0bbd/vvSlL43SiImIPty11167x2fW3LlzB7+ey+VwwQUXoK6uDhUVFTj11FOxfft2z88py2D6q1/9CpdeeimuueYavPLKK1i4cCGOO+44dHZ2jvbQiIhKeuqpp3DBBRfg+eefxxNPPAHbtvFP//RPSKfTu7U799xzsXXr1sF/brzxxlEaMRFRaQcccMBun1lPP/304NcuueQSPPzww3jooYfw1FNPoaOjA//6r//q+RllWS5qyZIlOPTQQ/Ef//EfAAAhBFpaWnDRRRfhqquuGuXRERF5s2PHDjQ2NuKpp57CUUcdBaAwY7po0SLcfvvtozs4IiIF1157LX7zm9/g1Vdf3eNriUQCDQ0N+MUvfoF/+7d/AwC88847mDdvHp577jl85CMfUX5O2c2YWpaFl19+Gccee+zga7qu49hjj8Vzzz03iiMjIhqaRCIBAKitrd3t9Z///Oeor6/HgQceiK9+9avIZDKjMTwiIiVr167F5MmTMWPGDJxxxhlob28HALz88suwbXu37DZ37lxMmzbNc3Yzh3XEw6Crqwuu62LSpEm7vT5p0iS88847ozQqIqKhEUJgxYoVOOKII3DggQcOvv6Zz3wGra2tmDx5Ml5//XVceeWVWL16NX7961+P4miJiD7YkiVLcN9992HOnDnYunUrrrvuOnz0ox/FqlWrsG3bNgSDQVRXV+/WZ9KkSdi2bZun55RdMCUiGk8uuOACrFq1are9WABw3nnnDf56/vz5aG5uxjHHHIP169dj5syZIz1MIqIPdcIJJwz+esGCBViyZAlaW1vx4IMPIhKJDNtzym4pv76+HoZh7HGSa/v27WhqahqlUREReXfhhRfikUcewcqVKzF16tQPbbtkyRIAwLp160ZiaERE+6S6uhr77bcf1q1bh6amJliWhb6+vt3aDCW7lV0wDQaDWLx4Mf785z8PviaEwJ///GcsXbp0FEdGRKRGSokLL7wQ//u//4snn3wSbW1tJfsUDxQ0Nzf7PDoion2XSqWwfv16NDc3Y/HixQgEArtlt9WrV6O9vd1zdivLpfxLL70UZ511Fg455BAcdthhuP3225FOp3H22WeP9tCIiEq64IIL8Itf/AK//e1vUVlZObjHqqqqCpFIBOvXr8cvfvELnHjiiairq8Prr7+OSy65BEcddRQWLFgwyqMnItrT5ZdfjpNOOgmtra3o6OjANddcA8MwcPrpp6OqqgrnnHMOLr30UtTW1iIej+Oiiy7C0qVLPZ3IB8o0mJ522mnYsWMHvvnNb2Lbtm1YtGgR/vCHP+xxIIqIqBzdeeedAAoloXb105/+FJ///OcRDAbxpz/9afCH7paWFpx66qn4+te/PgqjJSIqbfPmzTj99NPR3d2NhoYGHHnkkXj++efR0NAAALjtttug6zpOPfVU5PN5HHfccfjRj37k+TllWceUiIiIiCaesttjSkREREQTE4MpEREREZUFBlMiIiIiKgsMpkRERERUFhhMiYiIiKgsMJgSERERUVlgMCUiIiKissBgSkTko2uvvRaLFi0a7WF4Nn36dNx+++2jPQwimmAYTImIdvGXv/wFmqahr69vWN7v8ssv3+3+aCIi2jsGUyIiH0gp4TgOKioqUFdXt0/vZdv2sLYjIipXDKZENKYsX74cF154IS688EJUVVWhvr4e3/jGN7Dr7cq9vb0488wzUVNTg2g0ihNOOAFr164d/Pq7776Lk046CTU1NYjFYjjggAPw6KOPYtOmTTj66KMBADU1NdA0DZ///OcBAEII3HDDDWhra0MkEsHChQvx3//934PvWZxpfeyxx7B48WKEQiE8/fTTeyzlCyHwrW99C1OnTkUoFMKiRYvwhz/8YfDrmzZtgqZp+NWvfoVly5YhHA7j5z//+Qf+WWiahjvvvBMnn3wyYrEYrr/+eriui3POOWdwnHPmzMH3v//93fp9/vOfxymnnIKbb74Zzc3NqKurwwUXXPChwfY///M/UV1dzdlfIvKVOdoDICLy6v7778c555yDF198EX//+99x3nnnYdq0aTj33HMBFILX2rVr8bvf/Q7xeBxXXnklTjzxRLz11lsIBAK44IILYFkW/vrXvyIWi+Gtt95CRUUFWlpa8D//8z849dRTsXr1asTjcUQiEQDADTfcgAceeAB33XUXZs+ejb/+9a/47Gc/i4aGBixbtmxwbFdddRVuvvlmzJgxAzU1NfjLX/6y29i///3v45ZbbsGPf/xjHHTQQbj33ntx8skn480338Ts2bN3e59bbrkFBx10EMLh8F7/LK699lp873vfw+233w7TNCGEwNSpU/HQQw+hrq4Ozz77LM477zw0NzfjU5/61GC/lStXorm5GStXrsS6detw2mmnYdGiRYN/hru68cYbceONN+Lxxx/HYYcdNqT/z4iIlEgiojFk2bJlct68eVIIMfjalVdeKefNmyellHLNmjUSgHzmmWcGv97V1SUjkYh88MEHpZRSzp8/X1577bUf+P4rV66UAGRvb+/ga7lcTkajUfnss8/u1vacc86Rp59++m79fvOb3+zW5pprrpELFy4c/P3kyZPl9ddfv1ubQw89VJ5//vlSSik3btwoAcjbb7+95J8FALlixYqS7S644AJ56qmnDv7+rLPOkq2trdJxnMHXPvnJT8rTTjtt8Petra3ytttuk//n//wf2dzcLFetWlXyOURE+4ozpkQ05nzkIx+BpmmDv1+6dCluueUWuK6Lt99+G6ZpYsmSJYNfr6urw5w5c/D2228DAL7yla/gy1/+Mh5//HEce+yxOPXUU7FgwYK9Pm/dunXIZDL4+Mc/vtvrlmXhoIMO2u21Qw45ZK/v09/fj46ODhxxxBG7vX7EEUfgtddeU36fUu1++MMf4t5770V7ezuy2Swsy9qjMsABBxwAwzAGf9/c3Iw33nhjtza33HIL0uk0/v73v2PGjBlK4yEi2hfcY0pEE84Xv/hFbNiwAZ/73Ofwxhtv4JBDDsEdd9yx1/apVAoA8Pvf/x6vvvrq4D9vvfXWbvtMASAWiw3LGFXf5/3tfvnLX+Lyyy/HOeecg8cffxyvvvoqzj77bFiWtVu7QCCw2+81TYMQYrfXPvrRj8J1XTz44IND+DcgIvKOwZSIxpwXXnhht98///zzmD17NgzDwLx58+A4zm5turu7sXr1auy///6Dr7W0tOBLX/oSfv3rX+Oyyy7DT37yEwBAMBgEALiuO9h2//33RygUQnt7O2bNmrXbPy0tLcrjjsfjmDx5Mp555pndXn/mmWd2G9u+eOaZZ3D44Yfj/PPPx0EHHYRZs2Zh/fr1Q3qvww47DI899hi++93v4uabbx6W8RERfRgu5RPRmNPe3o5LL70U//7v/45XXnkFd9xxB2655RYAwOzZs/HP//zPOPfcc/HjH/8YlZWVuOqqqzBlyhT88z//MwBgxYoVOOGEE7Dffvuht7cXK1euxLx58wAAra2t0DQNjzzyCE488UREIhFUVlbi8ssvxyWXXAIhBI488kgkEgk888wziMfjOOuss5THfsUVV+Caa67BzJkzsWjRIvz0pz/Fq6++uteT917Nnj0b//Vf/4U//vGPaGtrw89+9jO89NJLaGtrG9L7HX744Xj00UdxwgknwDRNrFixYljGSUT0QRhMiWjMOfPMM5HNZnHYYYfBMAxcfPHFOO+88wa//tOf/hQXX3wxPvGJT8CyLBx11FF49NFHB5evXdfFBRdcgM2bNyMej+P444/HbbfdBgCYMmUKrrvuOlx11VU4++yzceaZZ+K+++7Dt7/9bTQ0NOCGG27Ahg0bUF1djYMPPhhf+9rXPI39K1/5ChKJBC677DJ0dnZi//33x+9+97vdTuTvi3//93/HP/7xD5x22mnQNA2nn346zj//fDz22GNDfs8jjzwSv//973HiiSfCMAxcdNFFwzJWIqL306TcpfgfEVGZW758ORYtWsTrMomIxiHuMSUiIiKissBgSkRERERlgUv5RERERFQWOGNKRERERGWBwZSIiIiIygKDKRERERGVBQZTIiIiIioLDKZEREREVBYYTImIiIioLDCYEhEREVFZYDAlIiIiorLAYEpEREREZeH/A+3lA2a4d2NiAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 800x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "f, ax = sbc_rank_plot(ranks, 1_000, plot_type=\"cdf\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above provides a visual representation of the cumulative density function (CDF) of `ranks` (blue and orange for each dimension of `theta`) with respect to the 95% confidence interval of a uniform distribution (grey)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# When things go haywire\n",
    "\n",
    "Next, we would like to explore some pathologies visible in sbc plots which can hint at our estimated posterior being somewhat wrong or completely off."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A shifted posterior mean\n",
    "\n",
    "In this scenario we emulate the situation that our posterior estimates incorrectly with a constant shift. We reuse our trained NPE posterior from above and wrap it so that all samples returned expose a constant shift by `+0.1`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "from utils_13_diagnosis_sbc import BiasedPosterior\n",
    "\n",
    "# this posterior shifts the expected value of the prior by .1\n",
    "posterior_ = BiasedPosterior(posterior, shift=0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6fe27c046ca14bd48780ae317f30d20c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Running 1000 sbc samples.:   0%|          | 0/1000 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'ks_pvals': tensor([0., 0.]), 'c2st_ranks': tensor([0.6815, 0.6730]), 'c2st_dap': tensor([0.5025, 0.4935])}\n"
     ]
    }
   ],
   "source": [
    "ranks, dap_samples = run_sbc(thetas, xs, posterior_)\n",
    "check_stats = check_sbc(ranks, thetas, dap_samples, 1_000)\n",
    "print(check_stats)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that the Kolmogorv-Smirnov p-values vanish (`'ks_pvals': tensor([0., 0.])`). Thus, we can reject the hypothesis that the `ranks` PDF is the uniform PDF. The `c2st` accuracies show values higher than `0.5`. This is indicative that the `ranks` distribution is not a uniform PDF as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdAAAAE9CAYAAAC7hzNcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAASvElEQVR4nO3de4xd1WEv4N/CGFyow8tuVGHKpGluSAGbh7mhNUGtXAFpIhyh2+g2bYmTVpVQH/emqiWkSCVCoaKt1UgNDZCECEiTVIBSk5KkwU3Ehb6oMaEmCQ5OgiOcUEhpwXmZazurf5xtMzYzZmbNzDln5nyfdDR79jl773UWXvz2WnvPXqXWGgBgeo4adAEAYD4SoADQQIACQAMBCgANBCgANBCgANDg6Ol8eNmyZXVsbGyOigILx9atW/+j1rp80OWYjLYMU3OktjytAB0bG8tDDz00O6WCBayU8s1Bl+FItGWYmiO1ZUO4ANBAgAJAAwEKAA2mdQ2U0bF3797s2rUre/bsGXRRhtqSJUuyYsWKLF68eNBFAfpMgDKhXbt2ZenSpRkbG0spZdDFGUq11jz77LPZtWtXXvWqVw26OECfGcJlQnv27Mkpp5wiPI+glJJTTjlFLx1GlABlUsLz5akjGF0ClHnjPe95TzZu3Dhr+9u5c2c+/vGPT/r+ZZddlhNPPDFvfvObZ+2YwMLhGihTcv+6dUd8/+K7757yvmqtqbXmqKMGe/52IEDf9ra3Tfj+hg0b8oMf/CA333xzn0sGzAd6oPTFzp0789rXvjZXXnllzjrrrDz55JO56qqrsnr16px55pm55pprDn52bGws11xzTc4777ycffbZ2b59+0v296EPfShvfOMb88Mf/vCQ9XfeeWfOOuusrFq1KhdffHGSZP/+/dmwYUMuuOCCrFy58mAgXn311XnggQdyzjnn5H3ve99LjrF27dosXbp0NqsBWED0QOmbHTt25LbbbsuFF16YJLnuuuty8sknZ//+/Vm7dm22bduWlStXJkmWLVuWhx9+OB/4wAeycePGfPjDHz64nxtuuCGbN2/Opk2bcuyxxx5yjGuvvTaf+9zncuqpp+a5555Lktxyyy054YQTsmXLlrzwwgtZs2ZNLrnkklx//fXZuHFj7rnnnv5UALCg9D1AJxoKnM7wH/PX6aeffjA8k+SOO+7IBz/4wezbty9PPfVUvvKVrxwM0CuuuCJJcv755+eTn/zkwW1uv/32nHbaadm0adOEf3u5Zs2arF+/Pm9961sP7uPee+/Ntm3bctdddyVJnn/++ezYsSPHHHPMnH3XUTDZsL72zKjQA6Vvjj/++IPLTzzxRDZu3JgtW7bkpJNOyvr16w/5c5ADPctFixZl3759B9efffbZeeSRRyb928ubbropDz74YD796U/n/PPPz9atW1Nrzfvf//5ceumlh3z2vvvum+VvCIwS10AZiN27d+f444/PCSeckKeffjqf/exnp7Tdueeem5tvvjmXX355vv3tb7/k/a9//et5/etfn2uvvTbLly/Pk08+mUsvvTQ33nhj9u7dmyR5/PHH8/3vfz9Lly7Nd7/73Vn9XsDo0ANlIFatWpVzzz03Z5xxRk477bSsWbNmyttedNFF2bhxY970pjdl8+bNWbZs2cH3NmzYkB07dqTWmrVr12bVqlVZuXJldu7cmfPOOy+11ixfvjybNm3KypUrs2jRoqxatSrr16/Pu971rkOO84Y3vCHbt2/P9773vaxYsSK33HLLS3qxwOgqtdYpf3j16tV1pnMIugY6Pzz22GN53eteN+hizAsT1VUpZWutdfWAivSy5qotJ9ozC8uR2rIhXABoIEABoIEABYAGApRJTef6+KhSRzC6BCgTWrJkSZ599lkBcQQH5gNdsmTJoIsCDIA/Y2FCK1asyK5du/Kd73xn0EUZakuWLMmKFSsGXQxgAAQoE1q8ePGET/oBoMcQLgA0EKAA0ECAAkADAQoADQQoADQQoADQQIACQAMBCgANBCgANBCgANBAgAJAAwEKAA0EKAA0EKAA0ECAAkADAQoADQQoADQQoADQQIACQAMBCgANBCgANBCgANBAgAJAAwEKAA0EKAA0EKAA0ECAAkADAQoADQQoADQQoADQQIACQAMBCgANBCgANBCgANBAgAJAAwEKAA0EKAA0EKAA0ECAAkADAQoADQQoADQQoADQQIACQAMBCgANBCgANBCgANBAgAJAAwEKAA0EKAA0EKAA0ECAAkADAQoADQQoADQQoADQQIACQAMBCgANBCgANBCgANBAgAJAAwEKAA0EKAA0EKAA0ECAAkADAQoADQQoADQQoADQQIACQAMBCgANBCgANBCgANBAgAJAAwEKAA0EKAA0EKAA0ECAAkADAQoADQQoADQQoADQQIACQAMBCgANBCgANBCgANBAgAJAAwEKAA0EKAA0EKAA0ECAAkADAQoADQQoADQQoADQQIACQAMBCgANBCgANBCgANBAgAJAAwEKAA0EKAA0EKAA0ECAAkADAQoADQQoADQQoADQQIACQAMBCgANBCgANBCgANBAgAJAAwEKAA0EKAA0EKAA0ECAAkADAQoADQQoADQ4etAFSJL7162bcP3Fd9/d55IAwNTogQJAAwEKAA0EKAA0EKAA0ECAAkADAQoADQQoADQQoADQQIACQAMBCgANBCgANBCgANBAgAJAAwEKAA2GYjozYOGYaHpCUxOyEOmBAkADAQoADQQoADQQoADQQIACQAMBCgANBCgANBCgANBAgAJAAwEKAA0EKAA0EKAA0ECAAkADAQoADQQoADQY6vlAJ5pXMDG3IACDpwcKAA0EKAA0EKAA0GCor4FOh+ulAPSTHigANBCgANBAgAJAAwEKAA0EKAA0EKAA0ECAAkADAQoADQQoADQQoADQQIACQIMF8yxcYGHwXGvmCz1QAGggQAGggQAFgAaugQJzznVNFqJ5GaCTNUYA6BdDuADQQIACQIN5OYQLLAwuxzCf6YECQAMBCgANBCgANBCgANBAgAJAAwEKAA0EKAA0EKAA0ECAAkADAQoADTzKD5gXTInGsNEDBYAGC74HOtFZqzNWAGZKDxQAGvS9B/rMiSdm76JF/T7sIbZs2TLQ4zPcjjvuuJx55pmDLgYw5PoeoHsXLcqx+/b1+7CHeMUrXjHQ4zPcdu/ePegiAPOAIVwAaCBAAaDBgr8LdyJfvu66Cdef+e5397kkMH8Nw/0MSfKJd77zJet+5qqrBlAShs1c388wkgEKzNww3M8wGfc5kMz9/QyGcAGggR7oFEw05Gu4F2C06YECQAMBCgANBCgANHANdJb5ExmA0aAHCgANBCgANBCgANBAgAJAAwEKAA3chQssOO6Gpx/0QAGggR4oMDIm65lORG+Vl6MHCgAN9ECHkOs3AMNPgAJMwDSGvBwBOs50ro8AtDDCtHC4BgoADfRA+8RwEMDCogcKAA30QOcRvViA4aEHCgAN9EAHyF2/ML+4g5bx9EABoIEABYAGAhQAGrgG2sj1S4DRJkAXKH/yAvOLNjv/GMIFgAYCFAAaGMKd56ZzLdbfsMHc0LZGU98DdPH+/XnhaLk9THbv3j3oIgyV4447btBFAOaBvifZTzz3XL8Pycu44IILBl0EgHlHVxBoYjTp5U00ujOdOjM6NDNzPZrkXz/QxGjSy5todOeH733vjLZneAhQgHnm/nXrJlx/8d1397kko02AAoygiUJYAE+PAGVSGhjA5AQowByZbKi1X9vP1vFm48R5IZ6QexIRADTQAwXgiBZi73E2CFDmlIYHLFSGcAGggR4o0zKXNzX42zaYmX7fdDTq9EABoIEeKM5agVkzSv8/EaAMvek0SMO9QL8IUAAGYr6fHAtQAKZtlIZqJyNA6bt+38k7jGeuwPznLlwAaCBAAaCBIVwA5q1BXrbRAwWABnqgAAy9YbzrV4Cy4E234blrF5gKQ7gA0KDUWqf+4VK+k+SbMzzmsiT/McN98CL1Oftmo05Pr7Uun43CzAVteWip09k1p215WgE6G0opD9VaV/f1oAuY+px96nRq1NPsU6eza67r0xAuADQQoADQYBAB+sEBHHMhU5+zT51OjXqafep0ds1pffb9GigALASGcAGgQd8CtJRyWSnlq6WUr5VSru7XcReCUsrOUsqjpZRHSikPdetOLqVsLqXs6H6e1K0vpZS/6Op5WynlvMGWfvBKKR8ppTxTSvnSuHXTrr9Sytu7z+8opbx9EN9lWGjPbbTlmRuq9lxrnfNXkkVJvp7kp5Mck+TfkvxsP469EF5JdiZZdti6P01ydbd8dZI/6ZZ/Oclnk5QkFyZ5cNDlH/QrycVJzkvypdb6S3Jykm90P0/qlk8a9HcbUH1qz+11py3PvA6Hpj33qwf6P5N8rdb6jVrr/0/y10mG78GG88u6JLd1y7clecu49bfXnn9JcmIp5ScHUL6hUWu9P8l/HrZ6uvV3aZLNtdb/rLX+V5LNSS6b88IPJ+15dmnL0zBM7blfAXpqkifH/b6rW8fU1CT3llK2llJ+u1v3ylrrU93yvyd5ZbesrqdmuvWnXl+kLtppy3NjIO3Zw+Tnh4tqrd8qpfxEks2llO3j36y11lKK26kbqT/6SFueY/2sw371QL+V5LRxv6/o1jEFtdZvdT+fSfI36Q2hPX1gOKf7+Uz3cXU9NdOtP/X6InXRSFueMwNpz/0K0C1JXlNKeVUp5Zgk/zvJp/p07HmtlHJ8KWXpgeUklyT5Unr1d+DOsbcnOTAH16eSXNndfXZhkufHDW3wounW3+eSXFJKOam7w++Sbt0o0p4baMtzajDtuY93Tv1yksfTu3vv3YO+k2u+vNK70/HfuteXD9RdklOSfD7JjiR/n+Tkbn1J8pddPT+aZPWgv8OgX0k+keSpJHvTu9bxmy31l+SdSb7Wvd4x6O814DrVnqdfZ9ry7NTj0LRnTyICgAaeRAQADQQoADQQoADQQIACQAMBCgANBOgMlFLeUkr52YbtLh/0DBallPtKKaun8flfKKXc0y3PqPyllN/tZkeopZRlrfuB2aQ9Nx/7Y93MPF/qZkpZ3Lqv+UaAzsxbkkyrwZVSjq61fqrWev10tml5b65Mt/wT+Mckv5Tkm7NUJJgNb4n23OJjSc5IcnaSH0vyW7NSsHlgJAO0lDJWStnenTk9Vkq5q5RyXPfe2lLKF0tvzr6PlFKO7dZfX0r5Sjen3MZSys8nuTzJn5Xe3H6v7l5/1z0o+oFSyhndtreWUm4qpTyY5E9LKetLKTeMK8sXuv1+vpTyUxNtc1j515dSPlVK+UKSz5dSfrzb9uGu3OvG7fuxUsqHSilfLqXcW0r5scP2dVR3rPdOUE+XdfX0cJIrDjv+gfLfWkq5sZTyL6WUb3Rnth/pjnvrRPVfa/1irXXn9P/LwUtpz4fsaxDt+TO1k+Rf03ss3mgY9FMlBvQki7H0ZkVY0/3+kSR/mGRJek/o/x/d+tuT/N/0nnLx1eTggydO7H7emuR/jdvv55O8plt+fZIvjPvcPUkWdb+vT3JDt/y3Sd4+7skYmyba5rDyr0/vCRwHnrZxdJJXdMvL0nuyRum+574k53Tv3ZHk17vl+9KbH+8TmeBJMuPq4jXdvu5Ics8E5b81vemsSnpTB+1O70z0qCRbDxx7kv8OO3PY3IheXtN9ac9D054XJ3k4yRsG/W+iX6+R7IF2nqy1/mO3/FdJLkry2iRP1Fof79bflt7krc8n2ZPkllLKFUl+cPjOSik/nuTnk9xZSnkkyc1Jxs/dd2etdf8E5fi5JB/vlj/alePltkm6uewOHD7JH5dStqX3GKtT8+J0Pk/UWh/plrem1wgPuDm9SWmvm2D/Z3Tb7qi91vFXk5QjSf62+8yjSZ6utT5aa/1Reo8rGzvCdjBbtOfBt+cPJLm/1vrAET6zoIxygB7+DMNJn2lYa92X3qwJdyV5c5K/m+BjRyV5rtZ6zrjX68a9//2GMh5pm/Hv/VqS5UnOr7Wek+Tp9M44k+SFcZ/bn0OnsPunJL9YSlmSmTlwjB8ddrwfxZR59If2PMD2XEq5pivzH8zw2PPKKAfoT5VSfq5bfluSf0hvWGeslPIz3frfSPL/urPRE2qtn0nyriSruve/m2RpktRadyd5opTyK0lSeg587kj+Kb3ZLJJew2k5ezshyTO11r2llF9McvoUt7slyWeS3FFeevPC9vTq4tXd77/aUC7oF+15QO25lPJbSS5N8qtdT3VkjHKAfjXJ75RSHktyUpIba617krwjvWGbR9M747opvUZ1Tzek8g958Szrr5Ns6G5SeHV6DeY3SykHZltYN4Vy/F6Sd3T7/o0k/6fhu3wsyequzFem11impNb650m+mOSjpZSjxq3fk+S3k3y6u+ngmUl2MW2llN8vpexK72aDbaWUD8/WvhlZ2nMG057Tq9NXJvnn7gasP5rFfQ+1kZyNpZQylt4F9LMGXRZgZrRnBmWUe6AA0Gwke6AAMFN6oADQQIACQAMBCgANBCgANBCgANBAgAJAg/8GGx0D2L5zTdcAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 576x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "f, ax = sbc_rank_plot(ranks, 1_000, plot_type=\"hist\", num_bins=30)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Inspecting the histograms for both dimenions, the rank distribution is clearly tilted to low rank values for both dimensions. Because we have shifted the expected value of the posterior to higher values (by `0.1`), we see more entries at low rank values."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's try to shift all posterior samples in the opposite direction. We shift the expectation value by `-0.1`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "posterior_ = BiasedPosterior(posterior, shift=-0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6a03e30496404b818bbb0c19cf1d6cc7",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Running 1000 sbc samples.:   0%|          | 0/1000 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'ks_pvals': tensor([0., 0.]), 'c2st_ranks': tensor([0.6795, 0.6955]), 'c2st_dap': tensor([0.4910, 0.4955])}\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdAAAAE9CAYAAAC7hzNcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAS1ElEQVR4nO3dfYxd5X0n8O/DS3FwHF5sN1phymSTLKQGmxezoTFBrVwBaVIcod1oy26JaatKUd+rtYQUqUQoVLS1GqmhAZIQAS1pBShrUpI0eBN1ocmGBRNq4uBgEhzhhAKlG2hwTDF5+sc9Y8Zmxp15fOfemTufj3Q1Z86995znPp6fv+c559xzSq01AMDMHDHsBgDAfCRAAaCBAAWABgIUABoIUABoIEABoMFRM3nxsmXL6tjY2Cw1BUbH1q1b/6nWunzY7ZiKWobpOVQtzyhAx8bG8uCDD/anVTDCSinfHXYbDkUtw/QcqpbtwgWABgIUABoIUABoMKNjoCwcL7/8cnbv3p29e/cOuylz2qJFi7JixYocffTRw24KMGAClEnt3r07S5YsydjYWEopw27OnFRrzXPPPZfdu3fnTW9607CbAwyYXbhMau/evVm6dKnwPIRSSpYuXWqUDguUAGVKwvPfp49g4RKgzBsf+tCHsmnTpr4tb9euXfn0pz895fMXX3xxjj/++LznPe/p2zqB0eEYKNNy7/r1h3z+grvumvayaq2pteaII4a7/TYeoJdddtmkz2/cuDF79uzJjTfeOOCWAfOBESgDsWvXrpx66qm5/PLLc/rpp+fJJ5/MBz7wgaxZsyYrV67MVVddtf+1Y2Njueqqq3L22WfnjDPOyI4dO16zvE984hN517velR/96EcHzL/jjjty+umnZ/Xq1bnggguSJK+88ko2btyYc889N6tWrdofiFdeeWXuu+++nHnmmfnIRz7ymnWsW7cuS5Ys6Wc3ACPECJSB2blzZ2655Zacd955SZJrrrkmJ554Yl555ZWsW7cu27Zty6pVq5Iky5Yty0MPPZSPfexj2bRpUz75yU/uX851112XLVu2ZPPmzTnmmGMOWMfVV1+dL37xiznppJPygx/8IEly00035bjjjssDDzyQl156KWvXrs2FF16Ya6+9Nps2bcrdd989mA4ARooAZWBOOeWU/eGZJLfffns+/vGPZ9++fXnqqafyzW9+c3+AXnrppUmSc845J5/5zGf2v+fWW2/NySefnM2bN0/63cu1a9dmw4YNed/73rd/Gffcc0+2bduWO++8M0ny/PPPZ+fOndn37LPZ9+KL+ZfHHz9gGUve8pb+fnBg1k12mGkmh5ZaCFAGZvHixfunn3jiiWzatCkPPPBATjjhhGzYsOGAr4OMjyyPPPLI7Nu3b//8M844Iw8//PCU37284YYbcv/99+dzn/tczjnnnGzdujW11nz0ox/NRRdddMBrP3/bbf3+iMAC4hgoQ/HCCy9k8eLFOe644/L000/nC1/4wrTed9ZZZ+XGG2/MJZdcku9///uvef7b3/523v72t+fqq6/O8uXL8+STT+aiiy7K9ddfn5dffjlJ8thjj+XFF1/M6xcvzg9ffLGvnwtYOIxAGYrVq1fnrLPOymmnnZaTTz45a9eunfZ7zz///GzatCnvfve7s2XLlixbtmz/cxs3bszOnTtTa826deuyevXqrFq1Krt27crZZ5+dWmuWL1+ezZs35/RTT82RRx6Zd/ziL+aySy/Nb15xxQHreec735kdO3bkhz/8YVasWJGbbrrpNaNYYOEqtdZpv3jNmjXVPQQXhkcffTRve9vbht2MWXXwsc9xMz0GOllflVK21lrXNDdulqllRs1sHQM9VC3bhQsADQQoADQQoADQQIAypZkcH1+o9BEsXAKUSS1atCjPPfecgDiE8fuBLlq0aNhNAYbA11iY1IoVK7J79+48++yzw27KrNn7zDOTzl/UfV90OhYtWpQVK1b0q0nAPCJAmdTRRx896ZV+Rsm9V1456fyzZvnyX8BosAsXABoIUABoIEABoIEABYAGAhQAGghQAGggQAGggQAFgAYCFAAaCFAAaCBAAaCBAAWABgIUABoIUABoIEABoIEABYAGAhQAGghQAGggQAGggQAFgAYCFAAaCFAAaCBAAaCBAAWABgIUABoIUABoIEABoIEABYAGAhQAGghQAGggQAGggQAFgAYCFAAaCFAAaCBAAaCBAAWABgIUABoIUABoIEABoIEABYAGAhQAGghQAGggQAGggQAFgAYCFAAaCFAAaCBAAaCBAAWABgIUABoIUABoIEABoIEABYAGAhQAGghQAGggQAGggQAFgAYCFAAaCFAAaCBAAaCBAAWABgIUABoIUABoIEABoIEABYAGAhQAGghQAGggQAGggQAFgAYCFAAaCFAAaCBAAaCBAAWABgIUABoIUABoIEABoIEABYAGAhQAGghQAGggQAGggQAFgAYCFAAaCFAAaCBAAaCBAAWABgIUABoIUABoIEABoIEABYAGRw27ATDb7l2/fthNAEaQESgANBCgANBAgAJAAwEKAA0EKAA0EKAA0ECAAkADAQoADQQoADQQoADQQIACQAMBCgANBCgANBCgANBAgAJAAwEKAA0EKAA0EKAA0ECAAkADAQoADQQoADQQoADQQIACQAMBCgANBCgANBCgANBAgAJAAwEKAA0EKAA0EKAA0ECAAkADAQoADQQoADQQoADQQIACQAMBCgANBCgANBCgANBAgAJAAwEKAA0EKAA0EKAA0OCoYTcAAKbr3vXrh92E/YxAAaCBAAWABgIUABoIUABoIEABoIEABYAGAhQAGghQAGggQAGggSsRMVLm0lVKgNFmBAoADYxAAZiT5voeJSNQAGggQAGggQAFgAYCFAAaCFAAaCBAAaCBAAWABgIUABoIUABoIEABoIEABYAGAhQAGghQAGggQAGggQAFgAYCFAAaCFAAaCBAAaCBAAWABgIUABoIUABoIEABoIEABYAGAhQAGghQAGhw1LAbAC3uXb9+2E0A+mS+1rMAZc6br8UFjDYBCsDAjNIGsWOgANBAgAJAA7twmTNGadcOLCST1e4Fd901hJYMlhEoADQQoADQQIACQAMBCgANBCgANBCgANDA11gA6LuF8LU0I1AAaGAEysAthC1TYPQZgQJAg4GPQLdv3549e/YMerXMId9bunTYTTik7du3Z+XKlcNuBjDHDTxA9+zZkze84Q2DXi1zyDH79g27CYdkAw+YDsdAgSb2JjFuru5Vmu29SQIUaGJvEuPm6l6l2d7AcxIRADQQoADQQIACQAMBCgANBCgANHAWLjOy/ZprJp2/8oMfHHBLAIbLCBQAGghQAGhgFy6zaqpdvgDznREoADQQoADQwC5cgBHjbPnBMAIFgAYCFAAa2IWL3T0ADYxAAaCBESjAAjaT72rbK3UgI1AAaCBAAaCBAAWABgIUABoIUABo4CxcgAXicO+O5O5KBxKg89xcuQiCwgIWGrtwAaCBEShTMqoEmJoRKAA0EKAA0ECAAkADx0DnkZkck5zstS4EDcMxV86Wp78EKMAc4w4p84MAnYOc/Qow9wlQgCGxsTy/DTxAjz322LzwwguDXu288tJRs/PPMlW/z9b65qtjjz122E0A5oGB/8+5cuXKQa9yzrp3/fpJ5580S+v70Yc/PND1zVf+RqfHxvD0zeZG6mT/BjaKe2Z7Y1gvA01saEzfVBuv/XDuuecOdH3zyWz/jfoeKAA0MAIF6KOpDs0weoxAAaCBEWijqbYyL7jrrhm9Hpg71CkzIUCBkTbTjd35RugPj124ANBgZEags7mVaQsPgIMZgQJAg5EZgQILx6gf12R+WJABqvgAOFwLMkABnNvA4ZqXAeoPH4Bhm5cBCoweG8bMNwK0z/wnAK8a9PkG6o9BGvkAVVAAzIaRD1BgMGysstC4kAIANJgTI9C58r1MW9AATNecCNCpCDQYTWqbUTDwAFU4AIwCx0ABoIEABYAGpdY6/ReX8myS7x7mOpcl+afDXAav0p/9148+PaXWurwfjZkNannO0qf9Nau1PKMA7YdSyoO11jUDXekI05/9p0+nRz/1nz7tr9nuT7twAaCBAAWABsMI0I8PYZ2jTH/2nz6dHv3Uf/q0v2a1Pwd+DBQARoFduADQYGABWkq5uJTyrVLK46WUKwe13lFQStlVSnmklPJwKeXBbt6JpZQtpZSd3c8TuvmllPJnXT9vK6WcPdzWD18p5VOllGdKKd+YMG/G/VdKeX/3+p2llPcP47PMFeq5jVo+fHOqnmuts/5IcmSSbyf5j0l+Isk/JPnpQax7FB5JdiVZdtC8P05yZTd9ZZI/6qZ/IckXkpQk5yW5f9jtH/YjyQVJzk7yjdb+S3Jiku90P0/opk8Y9mcbUn+q5/a+U8uH34dzpp4HNQL9z0ker7V+p9b6r0n+OomL4h6e9Ulu6aZvSfLeCfNvrT1fS3J8KeU/DKF9c0at9d4k/3zQ7Jn230VJttRa/7nW+v+TbEly8aw3fm5Sz/2llmdgLtXzoAL0pCRPTvh9dzeP6alJ7imlbC2l/Ho374211qe66X9M8sZuWl9Pz0z7T7++Sl+0U8uzYyj1PKdvZ8Z+59dav1dK+ckkW0opOyY+WWutpRSnUzfSfwyQWp5lg+zDQY1Av5fk5Am/r+jmMQ211u91P59J8r/S24X29PjunO7nM93L9fX0zLT/9Our9EUjtTxrhlLPgwrQB5K8tZTyplLKTyT5b0k+O6B1z2ullMWllCXj00kuTPKN9Ppv/Myx9ye5q5v+bJLLu7PPzkvy/IRdG7xqpv33xSQXllJO6M7wu7CbtxCp5wZqeVYNp54HeObULyR5LL2z9z447DO55ssjvTMd/6F7bB/vuyRLk3wpyc4k/zvJid38kuTPu35+JMmaYX+GYT+S/FWSp5K8nN6xjl9t6b8kv5Lk8e5xxbA/15D7VD3PvM/Ucn/6cc7UsysRAUADVyICgAYCFAAaCFAAaCBAAaCBAAWABgL0MJRS3ltK+emG910y7DtYlFL+rpSyZgav/9lSyt3d9GG1v5Tym93dEWopZVnrcqCf1HPzum/r7szzje5OKUe3Lmu+EaCH571JZlRwpZSjaq2frbVeO5P3tDw3W2ba/kl8JcnPJ/lun5oE/fDeqOcWtyU5LckZSV6X5Nf60rB5YEEGaCllrJSyo9tyerSUcmcp5djuuXWllK+X3j37PlVKOaabf20p5ZvdPeU2lVLekeSSJH9Sevf2e3P3+NvuQtH3lVJO6957cynlhlLK/Un+uJSyoZRy3YS2fLlb7pdKKT812XsOav+GUspnSylfTvKlUsrru/c+1LV7/YRlP1pK+UQpZXsp5Z5SyusOWtYR3bo+PEk/Xdz100NJLj1o/ePtv7mUcn0p5WullO90W7af6tZ782T9X2v9eq1118z/5eC11PMByxpGPX++dpL8v/Qui7cwDPuqEkO6ksVYendFWNv9/qkk/zPJovSu0P+fuvm3Jvnd9K5y8a1k/4Unju9+3pzkv0xY7peSvLWbfnuSL0943d1Jjux+35Dkum76b5K8f8KVMTZP9p6D2r8hvStwjF9t46gkb+iml6V3ZY3Sfc59Sc7snrs9yf/opv8uvfvj/VUmuZLMhL54a7es25PcPUn7b07vdlYlvVsHvZDelugRSbaOr3uKf4ddOejeiB4eM32o5zlTz0cneSjJO4f9NzGox4IcgXaerLV+pZv+yyTnJzk1yRO11se6+bekd/PW55PsTXJTKeXSJHsOXlgp5fVJ3pHkjlLKw0luTDLx3n131FpfmaQdP5Pk0930X3Tt+Pfek3T3shtffZI/LKVsS+8yVifl1dv5PFFrfbib3ppeEY67Mb2b0l4zyfJP6967s/aq4y+naEeS/E33mkeSPF1rfaTW+uP0Llc2doj3Qb+o5+HX88eS3Ftrve8QrxkpCzlAD76G4ZTXNKy17kvvrgl3JnlPkr+d5GVHJPlBrfXMCY+3TXj+xYY2Huo9E5/770mWJzmn1npmkqfT2+JMkpcmvO6VHHgLu68m+blSyqIcnvF1/Pig9f04bpnHYKjnIdZzKeWqrs2/f5jrnlcWcoD+VCnlZ7rpy5L8fXq7dcZKKW/p5v9ykv/TbY0eV2v9fJLfS7K6e/5fkixJklrrC0meKKX81yQpPeOvO5Svpnc3i6RXOC1bb8cleabW+nIp5eeSnDLN992U5PNJbi+vPXlhR3p98ebu919qaBcMinoeUj2XUn4tyUVJfqkbqS4YCzlAv5XkN0opjyY5Icn1tda9Sa5Ib7fNI+ltcd2QXlHd3e1S+fu8upX110k2dicpvDm9gvnVUsr43RbWT6Mdv5Xkim7Zv5zkdxo+y21J1nRtvjy9YpmWWuufJvl6kr8opRwxYf7eJL+e5HPdSQfPTLGIGSul/HYpZXd6JxtsK6V8sl/LZsFSzxlOPafXp29M8n+7E7D+oI/LntMW5N1YSilj6R1AP33YbQEOj3pmWBbyCBQAmi3IESgAHC4jUABoIEABoIEABYAGAhQAGghQAGggQAGgwb8By00B1HlyxwwAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 576x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ranks, dap_samples = run_sbc(thetas, xs, posterior_)\n",
    "check_stats = check_sbc(ranks, thetas, dap_samples, 1_000)\n",
    "print(check_stats)\n",
    "f, ax = sbc_rank_plot(ranks, 1_000, plot_type=\"hist\", num_bins=30)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A similar behavior is observed, but this time we see an upshot of ranks to higher values. Because we have shifted the expected value of the posterior to smaller values, we see an upshot in high rank counts.\n",
    "\n",
    "It is interesting to see that the historgams obtained provide very convincing evidence that this is not a uniform distribution. \n",
    "\n",
    "To conlude at this point, **the rank distribution is capable of identifying pathologies of the estimated posterior**:\n",
    "\n",
    "+ a **left-skewed rank distribution** shows a systematic **underestimation of the posterior mean**  \n",
    "(we shifted the posterior by `0.1`)\n",
    "+ a **rank-skewed rank distribution** shows a systematic **overestimation of the posterior mean**  \n",
    "(we shifted the posterior by `-0.1`)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A dispersed posterior\n",
    "\n",
    "In this scenario we emulate the situation if our posterior estimates incorrectly with a dispersion, i.e. the posterior is too wide or too thin. We reuse our trained NPE posterior from above and wrap it so that all samples return a dispersion by 100% more wide (`2`), i.e. the variance is overestimated by a factor of 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "from utils_13_diagnosis_sbc import DispersedPosterior\n",
    "\n",
    "# this posterior which disperses the expected posterior value of the prior by 2.\n",
    "posterior_ = DispersedPosterior(posterior, dispersion=2.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "5fbe9fcd9c5e4e4ea0d3fa90149cc756",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Running 1000 sbc samples.:   0%|          | 0/1000 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'ks_pvals': tensor([8.2151e-09, 6.9635e-07]), 'c2st_ranks': tensor([0.6150, 0.6160]), 'c2st_dap': tensor([0.5050, 0.4905])}\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdAAAAE9CAYAAAC7hzNcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAATI0lEQVR4nO3df4xd1WEn8O/hp2vq8MP2RitMcTZNQ2qw+WE2tCZoK1dAmghHqI3atCVOWlVC2+5uVmsJydISobAiXWsjbdgASYiANkkFKGunJGlwE1Fou2WNCTUhODgJjnDCQko3OIljFpOzf7xrM7ZnzLwzb+bNvPl8pKe5c9+79553PMffd+4979xSaw0A0J/jhl0AAJiLBCgANBCgANBAgAJAAwEKAA0EKAA0OKGfFy9ZsqQuX758mooCo2P79u3/VGtdOuxyTERbhsk5VlvuK0CXL1+eRx55ZDClghFWSvnusMtwLNoyTM6x2rJTuADQQIACQAMBCgAN+roGyvzx8ssvZ8+ePdm/f/+wizKrLViwIMuWLcuJJ5447KIAM0yAMq49e/Zk0aJFWb58eUopwy7OrFRrzQsvvJA9e/bkDW94w7CLA8wwp3AZ1/79+7N48WLheQyllCxevFgvHeYpAcqEhOdrU0cwfwlQ5owPfvCD2bRp08D2t3v37nzmM5+Z8Pkrr7wyp512Wt75zncO7JjA6HANlEl5cN26Yz5/2ZYtk95XrTW11hx33HA/vx0M0Pe85z3jPr9hw4bs27cvt9122wyXDJgL9ECZEbt3786b3/zmXHPNNTn33HPzzDPP5Nprr83q1auzYsWKXH/99Ydeu3z58lx//fW58MILc95552Xnzp1H7e8Tn/hE3v72t+enP/3pYevvueeenHvuuVm1alUuu+yyJMkrr7ySDRs25OKLL87KlSsPBeJ1112Xhx56KOeff34+8pGPHHWMtWvXZtGiRYOsBmCE6IEyY3bt2pU777wzl1xySZLkxhtvzBlnnJFXXnkla9euzY4dO7Jy5cokyZIlS/Loo4/mYx/7WDZt2pRPfvKTh/Zz8803Z+vWrdm8eXNOPvnkw45xww035Mtf/nLOPPPM/PCHP0yS3H777Tn11FOzbdu2vPTSS1mzZk0uv/zy3HTTTdm0aVPuu+++makAYKQI0DliolOo/Zw6Hbazzz77UHgmyd13352Pf/zjOXDgQJ599tl84xvfOBSgV199dZLkoosuyuc+97lD29x1110566yzsnnz5nG/e7lmzZqsX78+7373uw/t4/7778+OHTty7733JklefPHF7Nq1KyeddNK0vVdoMQrtfD4RoMyYU0455dDy008/nU2bNmXbtm05/fTTs379+sO+DnKwZ3n88cfnwIEDh9afd955eeyxxyb87uWtt96ahx9+OF/4whdy0UUXZfv27am15qMf/WiuuOKKw177wAMPDPgdAvOJa6AMxd69e3PKKafk1FNPzXPPPZcvfelLk9ruggsuyG233Zarrroq3//+9496/tvf/nbe+ta35oYbbsjSpUvzzDPP5Iorrsgtt9ySl19+OUny1FNP5Sc/+UkWLVqUH/3oRwN9X8D8oQfKUKxatSoXXHBBzjnnnJx11llZs2bNpLe99NJLs2nTprzjHe/I1q1bs2TJkkPPbdiwIbt27UqtNWvXrs2qVauycuXK7N69OxdeeGFqrVm6dGk2b96clStX5vjjj8+qVauyfv36fOADHzjsOG9729uyc+fO/PjHP86yZcty++23H9WLBeavUmud9ItXr15d3UNwOGb62siTTz6Zt7zlLdOy71EzXl2VUrbXWlcPqUivSVuenVwDnX2O1Zb1QOc4DQ5gOFwDBYAGAhQAGghQJtTP9fH5Sh3B/CVAGdeCBQvywgsvCIhjOHg/0AULFgy7KMAQGETEuJYtW5Y9e/bkBz/4wbCLMqstWLAgy5YtG3YxgCEQoIzrxBNPHHemHwB6nMIFgAYCFAAaOIULMAeZRGX49EABoIEABYAGAhQAGghQAGhgENGAubAPMD/ogQJAAwEKAA0EKAA0EKAA0MAgIoAhmGjAIXOHHigANBCgANBAgAJAA9dAZ6HZcm2kn3KYKAKYbwQowDSaLR+IGTyncAGggQAFgAYCFAAauAYKNHniiSeyb9++YRdj1vve4sVT3se2bdsmvd/xXjtfLVy4MCtWrJi2/QtQoMm+ffvyute9btjFmPVOPnBgyvsYr54n2q9/k1ft3bt3WvfvFC4ANBCgANBAgAJAgxm/BjrqAw8GcWF/ugYd9Kufcnz2/e8/at0vXnvtlMswDNM98AAYDTMeoKM+8GAQF/ana9BBv6Zajrn67zzdAw+A0eAULgA0EKAA0ECAAkADAQoADQQoADQQoADQQIACQAMBCgAN3I0FYECeuPHGYRdh3DKs2LhxCCUZfXqgANBAgAJAAwEKAA1cAx1R/VwHma7rNhPt1/UYYBQIUIARN4gPswYnHc0pXABoIEABoIFTuJPg1AUARxKgwMjy4Zfp5BQuADQQoADQQIACQAPXQIdoNkw8DUAbAQpwDAYiMRGncAGggQAFgAYCFAAaCFAAaCBAAaCBAAWABgIUABoIUABoYCKFMfqZGcgsQgDzmwAF5pWJPvz2M7uQD9AkTuECQBMBCgANBCgANBCgANBAgAJAAwEKAA0EKAA0EKAA0MBECsCcN+oTG4z6+5ur9EABoIEABYAGAhQAGrgGOo+4jgIwOHqgANBAgAJAAwEKAA0EKAA0EKAA0MAoXAAOMVp/8vRAAaCBAAWABk7hMuMmOkW0YuPGGS4JQDs9UABooAcKzCkGuQyOupwaPVAAaCBAAaCBAAWABgIUABoIUABoIEABoIEABYAGAhQAGphIAYAm/UzEMIpTdeqBAkADAQoADeblKVzzPwIwVXqgANBAgAJAAwEKAA0EKAA0GPlBRAYMATAd9EABoIEABYAGAhQAGghQAGggQAGggQAFgAYj/zUWgMnwlTf6pQcKAA0EKAA0EKAA0ECAAkADAQoADQQoADTwNZYZYoj8axuvjlZs3DiEkgC8Nj1QAGigBwrArDLRGbvZdkZKDxQAGghQAGggQAGggQAFgAYGEQFD5ytMzEV6oADQQA+UWW2uDGcH5h89UABoIEABoIFTuMCsZP7o0TKKl2P0QAGggQAFgAYCFAAaCFAAaCBAAaCBAAWABgIUABoIUABoIEABoIEABYAGAhQAGghQAGggQAGggQAFgAYCFAAaCFAAaCBAAaCBAAWABgIUABoIUABocMKwC3AsT9x447jrV2zcOOnXAsB00AMFgAYCFAAazOpTuMBocamFUaIHCgAN9ECZk/oZYAYwHfRAAaCBAAWABk7hAk0WLlyYvXv39rXNSyf4L4fDPfrhDx+9coK/k37/3hYuXNhSpEmb8b/mfhrdRI1tvO01TJL+G9h4prvRjYoVK1b0vc1PP/ShaSgJ88XFF1887CIcZsZTp59GN1FjG68SNUyS2dfAgNHlGigANBCgANDAhUMA5oQH1607at1lW7YMoSQ9czJAx6tESGZfAwNGl1O4ANBAgAJAAwEKAA0EKAA0EKAA0ECAAkADAQoADQQoADSYkxMpQD8mmnjDBAvAVOiBAkADAQoADQQoADQQoADQQIACQAMBCgANBCgANBCgANBAgAJAAwEKAA0EKAA0EKAA0MBk8sCUTDRZP4w6PVAAaCBAAaCBAAWABgIUABoIUABoIEABoIEABYAGAhQAGghQAGggQAGgwayYys9UYADMNXqgANBAgAJAAwEKAA0EKAA0EKAA0ECAAkADAQoADQQoADQQoADQQIACQAMBCgANBCgANBCgANBAgAJAAwEKAA0EKAA0EKAA0OCEYRcAAFo9uG7duOsv27Jl2o+tBwoADQQoADQQoADQQIACQAMBCgANBCgANBCgANBAgAJAAwEKAA1mfCaiiWaNAIC5RA8UABoIUABoIEABoIEABYAGAhQAGghQAGggQAGggQAFgAYCFAAaCFAAaCBAAaCBAAWABgIUABrM+N1YYLYY785Al23ZMoSSAHORHigANBCgANBAgAJAAwEKAA0EKAA0EKAA0ECAAkADAQoADQQoADQotdbJv7iUHyT57hSPuSTJP01xH7xKfQ7eIOr07Frr0kEUZjpoy7OWOh2saW3LfQXoIJRSHqm1rp7Rg44w9Tl46nRy1NPgqdPBmu76dAoXABoIUABoMIwA/fgQjjnK1OfgqdPJUU+Dp04Ha1rrc8avgQLAKHAKFwAazFiAllKuLKV8s5TyrVLKdTN13FFQStldSnm8lPJYKeWRbt0ZpZStpZRd3c/Tu/WllPLfu3reUUq5cLilH75SyqdKKc+XUr4+Zl3f9VdKeW/3+l2llPcO473MFtpzG2156mZVe661TvsjyfFJvp3kXyU5Kck/JvnlmTj2KDyS7E6y5Ih1f5rkum75uiQf7pZ/I8mXkpQklyR5eNjlH/YjyWVJLkzy9db6S3JGku90P0/vlk8f9nsbUn1qz+11py1PvQ5nTXueqR7ov07yrVrrd2qt/y/JXyRZN0PHHlXrktzZLd+Z5F1j1t9Ve/4hyWmllH85hPLNGrXWB5P88xGr+62/K5JsrbX+c631/ybZmuTKaS/87KQ9D5a23IfZ1J5nKkDPTPLMmN/3dOuYnJrk/lLK9lLKH3XrXl9rfbZb/j9JXt8tq+vJ6bf+1Our1EU7bXl6DKU9n9B/ORmCS2ut3yul/IskW0spO8c+WWutpRTDqRupP2aQtjzNZrIOZ6oH+r0kZ435fVm3jkmotX6v+/l8kv+Z3im05w6ezul+Pt+9XF1PTr/1p15fpS4aacvTZijteaYCdFuSN5VS3lBKOSnJbyf5/Awde04rpZxSSll0cDnJ5Um+nl79HRw59t4kW7rlzye5pht9dkmSF8ec2uBV/dbfl5NcXko5vRvhd3m3bj7Snhtoy9NqOO15BkdO/UaSp9Ibvbdx2CO55sojvZGO/9g9njhYd0kWJ/lKkl1J/jrJGd36kuR/dPX8eJLVw34Pw34k+WySZ5O8nN61jj9oqb8k70/yre7xvmG/ryHXqfbcf51py4Opx1nTns1EBAANzEQEAA0EKAA0EKAA0ECAAkADAQoADQToFJRS3lVK+eWG7a4a9h0sSikPlFJW9/H6f1NKua9bnlL5Syl/3N0doZZSlrTuBwZJe24+9qe7O/N8vbtTyomt+5prBOjUvCtJXw2ulHJCrfXztdab+tmm5bnp0m/5x/F3SX49yXcHVCQYhHdFe27x6STnJDkvyc8l+cOBFGwOmJcBWkpZXkrZ2X1yerKUcm8pZWH33NpSytdK7559nyqlnNytv6mU8o3unnKbSim/muSqJP+19O7t98bu8VfdRNEPlVLO6ba9o5Ryaynl4SR/WkpZX0q5eUxZvtrt9yullF8Yb5sjyr++lPL5UspXk3yllPLz3baPduVeN2bfT5ZSPlFKeaKUcn8p5eeO2Ndx3bE+NE49XdnV06NJrj7i+AfLf0cp5ZZSyj+UUr7TfbL9VHfcO8ar/1rr12qtu/v/l4Ojac+H7WsY7fmLtZPkf6c3Ld78MOxZJYY0k8Xy9O6KsKb7/VNJ/lOSBenN0P9L3fq7kvyH9Ga5+GZyaOKJ07qfdyT5zTH7/UqSN3XLb03y1TGvuy/J8d3v65Pc3C3/ZZL3jpkZY/N42xxR/vXpzcBxcLaNE5K8rltekt7MGqV7nweSnN89d3eS3+uWH0jv/nifzTgzyYypizd1+7o7yX3jlP+O9G5nVdK7ddDe9D6JHpdk+8FjT/DvsDtH3BvRw6Pfh/Y8a9rziUkeTfK2Yf9NzNRjXvZAO8/UWv+uW/7zJJcmeXOSp2utT3Xr70zv5q0vJtmf5PZSytVJ9h25s1LKzyf51ST3lFIeS3JbkrH37run1vrKOOX4lSSf6Zb/rCvHa22TdPeyO3j4JP+llLIjvWmszsyrt/N5utb6WLe8Pb1GeNBt6d2U9sZx9n9Ot+2u2msdfz5BOZLkL7vXPJ7kuVrr47XWn6U3XdnyY2wHg6I9D789fyzJg7XWh47xmpEynwP0yDkMJ5zTsNZ6IL27Jtyb5J1J/mqclx2X5Ie11vPHPN4y5vmfNJTxWNuMfe53kyxNclGt9fwkz6X3iTNJXhrzuldy+C3s/j7Jr5VSFmRqDh7jZ0cc72dxyzxmhvY8xPZcSrm+K/N/nOKx55T5HKC/UEr5lW75PUn+Nr3TOstLKb/Yrf/9JH/TfRo9tdb6xSQfSLKqe/5HSRYlSa11b5KnSym/lSSl5+DrjuXv07ubRdJrOC2f3k5N8nyt9eVSyq8lOXuS292e5ItJ7i5HD17YmV5dvLH7/XcaygUzRXseUnsupfxhkiuS/E7XU5035nOAfjPJvy2lPJnk9CS31Fr3J3lfeqdtHk/vE9et6TWq+7pTKn+bVz9l/UWSDd0ghTem12D+oJRy8G4L6yZRjj9J8r5u37+f5N83vJdPJ1ndlfma9BrLpNRa/1uSryX5s1LKcWPW70/yR0m+0A06eH6CXfStlPLvSil70htssKOU8slB7Zt5S3vOcNpzenX6+iT/qxuA9Z8HuO9ZbV7ejaWUsjy9C+jnDrsswNRozwzLfO6BAkCzedkDBYCp0gMFgAYCFAAaCFAAaCBAAaCBAAWABgIUABr8f2qKIehQZM5RAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 576x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ranks, dap_samples = run_sbc(thetas, xs, posterior_)\n",
    "check_stats = check_sbc(ranks, thetas, dap_samples, 1_000)\n",
    "print(check_stats)\n",
    "f, ax = sbc_rank_plot(ranks, 1_000, plot_type=\"hist\", num_bins=30)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The rank histograms now look more like a very wide gaussian distribution centered in the middle. The KS p-values again vanish unsurprisingly (we must reject the hypothesis that both distributions are from the same uniform PDF) and the c2st_ranks indicate that the rank histogram is not uniform too. As our posterior samples are distributed too broad now, we obtain more \"medium\" range ranks and hence produce the peak of ranks in the center of the histogram."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can repeat this exercise by making our posterior too thin, i.e. the variance of the posterior is too small. Let's cut it by half."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "posterior_ = DispersedPosterior(posterior, dispersion=0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ff39afff13e8495cab99ae88c5b588ee",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Running 1000 sbc samples.:   0%|          | 0/1000 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'ks_pvals': tensor([1.1753e-07, 1.7929e-08]), 'c2st_ranks': tensor([0.5755, 0.6125]), 'c2st_dap': tensor([0.4980, 0.5075])}\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdAAAAE9CAYAAAC7hzNcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAATJ0lEQVR4nO3df4xd1WEn8O/BuLhxDAbbjVaY4mySQmqw+WE2NCZoK1dAmhRHqI1atiWmrSpF7e62q1pCirREKKxI19pIGzZAEiKgJa0AZZ2UJA1uohbabqkxoSYODibBEU5YTGmBBses7Zz9412bsRk7M8dv3puZ9/lIT3PnvvvuPe94jr/3nPur1FoDAEzOCcMuAADMRAIUABoIUABoIEABoIEABYAGAhQAGpw4mYUXL15cly1bNkVFgdljy5Yt/1RrXTLschyNtgwTc6y2PKkAXbZsWR555JH+lApmsVLKd4ddhmPRlmFijtWWDeECQAMBCgANBCgANJjUMVBGx759+7Jr167s3bt32EWZ1ubNm5elS5dm7ty5wy4KMGAClHHt2rUrCxYsyLJly1JKGXZxpqVaa1544YXs2rUrb37zm4ddHGDADOEyrr1792bRokXC8xhKKVm0aJFeOowoAcpRCc8fTx3B6BKgzBgf/vCHs2HDhr6tb+fOnfnsZz971PevuOKKLFy4MO9973v7tk1g9nAMlAl5cO3aY75/6ec/P+F11VpTa80JJwx3/+1ggF599dXjvr9+/frs2bMnt91224BLBswEeqAMxM6dO3PWWWflmmuuyTnnnJNnnnkmH/zgB7Nq1aosX748119//aFlly1bluuvvz4XXHBBzj333Gzfvv116/vUpz6Vd7/73fnhD3942Px7770355xzTlauXJlLL700SXLgwIGsX78+F110UVasWHEoEK+77ro89NBDOe+88/Kxj33sddtYs2ZNFixY0M9qAGYRPVAGZseOHbnzzjtz8cUXJ0luvPHGnHbaaTlw4EDWrFmTrVu3ZsWKFUmSxYsX59FHH80nPvGJbNiwIZ/+9KcPrefmm2/Opk2bsnHjxpx00kmHbeOGG27IV77ylZx++ul58cUXkyS33357TjnllGzevDmvvvpqVq9encsuuyw33XRTNmzYkPvvv38wFQDMKgKUgTnzzDMPhWeS3HPPPfnkJz+Z/fv359lnn803v/nNQwF61VVXJUkuvPDCfO5znzv0mbvuuitnnHFGNm7cOO61l6tXr866devy/ve//9A6HnjggWzdujX33XdfkuSll17Kjh07sv/557P/lVfyr089ddg6Frz1rf394sCUG+8w02QOLbUQoAzM/PnzD00//fTT2bBhQzZv3pxTTz0169atO+xykIM9yzlz5mT//v2H5p977rl57LHHjnrt5a233pqHH344X/ziF3PhhRdmy5YtqbXm4x//eC6//PLDlv3S3Xf3+ysCI8QxUIbi5Zdfzvz583PKKafkueeey5e//OUJfe7888/PbbfdliuvvDLf//73X/f+t7/97bzjHe/IDTfckCVLluSZZ57J5ZdfnltuuSX79u1Lkjz55JN55ZVX8sb58/ODV17p6/cCRsfAe6DD6GYz/axcuTLnn39+zj777JxxxhlZvXr1hD97ySWXZMOGDXnPe96TTZs2ZfHixYfeW79+fXbs2JFaa9asWZOVK1dmxYoV2blzZy644ILUWrNkyZJs3Lgx55x1VubMmZN3/tIv5eqrrsrvXXvtYdt517vele3bt+cHP/hBli5dmttvv/11vVhgdJVa64QXXrVqVT3eZwgK0JnhiSeeyNvf/vZhF2NKHXns86DJHgMdr65KKVtqrauaCzfF+tGWYTqZqmw5Vls2hAsADQQoADQQoADQQIByVJM5Pj6q1BGMLgHKuObNm5cXXnhBQBzDweeBzps3b9hFAYbAjRQY19KlS7Nr1648//zzwy7KlNm7e/e48+d114tOxLx587J06dJ+FQmYQQQo45o7d+64d/qZTR687rpx55/vsipgAgQo0ORoj7hzXTejwjFQAGggQAGggQAFgAYCFAAaCFAAaCBAAaCBAAWABgIUABoIUABoIEABoIEABYAGAhQAGghQAGggQAGggQAFgAYCFAAaCFAAaCBAAaCBAAWABgIUABoIUABoIEABoIEABYAGAhQAGghQAGggQAGggQAFgAYCFAAaCFAAaCBAAaCBAAWABgIUABoIUABoIEABoIEABYAGAhQAGghQAGggQAGggQAFgAYCFAAaCFAAaCBAAaCBAAWABgIUABoIUABoIEABoIEABYAGAhQAGghQAGggQAGggQAFgAYCFAAaCFAAaCBAAaCBAAWABgIUABoIUABoIEABoIEABYAGAhQAGghQAGggQAGggQAFgAYCFAAaCFAAaCBAAaCBAAWABgIUABoIUABoIEABoIEABYAGAhQAGghQAGggQAGggQAFgAYCFAAaCFAAaCBAAaCBAAWABgIUABoIUABoIEABoIEABYAGAhQAGghQAGggQAGggQAFgAYCFAAaCFAAaCBAAaCBAAWABgIUABoIUABoIEABoIEABYAGAhQAGghQAGggQAGggQAFgAYCFAAaCFAAaCBAAaCBAAWABgIUABoIUABoIEABoIEABYAGAhQAGghQAGggQAGggQAFgAYCFAAanDjoDe5euDD75sw5bN7mzZsHXQxGyFO33DL+G4sWjTt727ZtWb58+RSWCJgNBh6g++bMyUn79x827+STTx50MRghR/69/Th79uyZopIAs4khXABoIEABoIEABYAGAhQAGghQAGggQAGggQAFgAYDvw4UptK2G28cdhGAEaEHCgANBCgANBCgANDAMVAAZowH164ddhEOmdYBerQTQpZ/6EMDLgkAHM4QLgA0EKAA0GBaDOG6dg+AmWZaBCgw8+xeuDD75sx53fzNmzcPoTSMiu8tWjThZbdt25bly5dPWVkEKNBk35w5OWn//tfNP/nkk4dQGkbFeH9zR7Nnz54pLMmIBqizewE4Xk4iAoAGAhQAGozkEC7Tk6F1YCbRAwWABgIUABoIUABoIEABoIGTiIAp5wQxZiM9UABoIEABoIEABYAGAhQAGjiJCOgrz/dlVOiBAkADAQoADQQoADSYkcdAxzvGMpsuyO7HReezvY4Ahm1GBigAs9tMOBnNEC4ANNADneFmwl7a8RqF7ziq3COXmUwPFAAaCFAAaDDrh3AN/wEwFfRAAaCBAAWABrN+CJfpx7A6MBvogQJAAz1Q+sL1fMCo0QMFgAZ6oADThIdAzCx6oADQQIACQAMBCgANZs0xUNcWAjBIeqAA0ECAAkCDWTOEC8xu0/1mHS5B+fGm+7/hZOmBAkADAQoADQzhMinOdobRNNuGX/tBDxQAGuiBAjBUM3VkS4AOyEz9AwFgfIZwAaCBHigw7fRjxMZ1mYMxyvWsBwoADfRAATjE+RoTpwcKAA0EKAA0EKAA0MAx0AkY1bPMHAsBODoBCjCO2bQD6T62U8MQLgA00ANlSs2mvXiAsQQowDQ2lcOvU7WDOyo7zoZwAaCBHmgjB+WBVqPSQ5vt9EABoIEABYAGAx/CnXvgQF49cXqOHL/88svjzp9MeR/96EfHf6MP33m88k2mbEf7fkczXf+dptob3vCGYRcBmAEG/j/kT7344qA3OWEXXXTRuPN/+JGPDLgk4xuvfJMp29G+39FMl+89aMuXLx92EYAZYDS7GMBIcvIf/SRAyYNr1w67CAAzjgAFmCJTebmKS2GGT4COED1NgP5xGQsANNADBZpMl0vS+nHp2PFeIsb0NNWXpPkLAZpM50vSJut4LxFjeprqS9IEKDDynB9AC8dAAaCBHugMYi8ZYPrQAwWABgIUABoYwh3DECkAE6UHCgANBCgANBCgANBAgAJAAwEKAA0EKAA0EKAA0ECAAkADAQoADQQoADQQoADQQIACQAMBCgANBCgANBCgANBAgAJAAwEKAA0EKAA0EKAA0ECAAkADAQoADQQoADQQoADQQIACQAMBCgANBCgANBCgANBAgAJAAwEKAA0EKAA0EKAA0ECAAkADAQoADQQoADQQoADQQIACQAMBCgANBCgANBCgANCg1FonvnApzyf57nFuc3GSfzrOdfAa9dl//ajTM2utS/pRmKmgLU9b6rS/prQtTypA+6GU8kitddVANzqLqc/+U6cTo576T53211TXpyFcAGggQAGgwTAC9JND2OZspj77T51OjHrqP3XaX1NanwM/BgoAs4EhXABoMLAALaVcUUr5VinlqVLKdYPa7mxQStlZSnm8lPJYKeWRbt5ppZRNpZQd3c9Tu/mllPI/u3reWkq5YLilH75SymdKKbtLKd8YM2/S9VdK+UC3/I5SygeG8V2mC+25jbZ8/KZVe661TvkryZwk307yb5P8RJJ/TPKzg9j2bHgl2Zlk8RHz/ijJdd30dUk+2k3/YpIvJylJLk7y8LDLP+xXkkuTXJDkG631l+S0JN/pfp7aTZ867O82pPrUntvrTls+/jqcNu15UD3Qf5fkqVrrd2qt/y/JnyVZO6Btz1Zrk9zZTd+Z5H1j5t9Ve/4+ycJSyr8ZQvmmjVrrg0n++YjZk62/y5NsqrX+c631X5JsSnLFlBd+etKe+0tbnoTp1J4HFaCnJ3lmzO+7unlMTE3yQCllSynld7p5b6q1PttN/98kb+qm1fXETLb+1Otr1EU7bXlqDKU9nzj5cjIEl9Rav1dK+akkm0op28e+WWutpRSnUzdSfwyQtjzFBlmHg+qBfi/JGWN+X9rNYwJqrd/rfu5O8r/TG0J77uBwTvdzd7e4up6Yydafen2NumikLU+ZobTnQQXo5iRvK6W8uZTyE0l+NckXBrTtGa2UMr+UsuDgdJLLknwjvfo7eObYB5J8vpv+QpJrurPPLk7y0pihDV4z2fr7SpLLSimndmf4XdbNG0XacwNteUoNpz0P8MypX0zyZHpn731o2GdyzZRXemc6/mP32naw7pIsSvLVJDuS/GWS07r5Jcn/6ur58SSrhv0dhv1K8qdJnk2yL71jHb/VUn9JfjPJU93r2mF/ryHXqfY8+TrTlvtTj9OmPbsTEQA0cCciAGggQAGggQAFgAYCFAAaCFAAaCBAj0Mp5X2llJ9t+NyVw36CRSnlr0opqyax/L8vpdzfTR9X+Uspv9c9HaGWUha3rgf6SXtu3vbd3ZN5vtE9KWVu67pmGgF6fN6XZFINrpRyYq31C7XWmybzmZb3pspkyz+Ov03yC0m+26ciQT+8L9pzi7uTnJ3k3CQ/meS3+1KwGWAkA7SUsqyUsr3bc3qilHJfKeUN3XtrSilfL71n9n2mlHJSN/+mUso3u2fKbSilvDPJlUn+e+k92+8t3esvuhtFP1RKObv77B2llFtLKQ8n+aNSyrpSys1jyvK1br1fLaX89HifOaL860opXyilfC3JV0spb+w++2hX7rVj1v1EKeVTpZRtpZQHSik/ecS6Tui29ZFx6umKrp4eTXLVEds/WP47Sim3lFL+vpTynW7P9jPddu8Yr/5rrV+vte6c/L8cvJ72fNi6htGev1Q7Sf4hvdvijYZh31ViSHeyWJbeUxFWd79/JskfJpmX3h36f6abf1eS30/vLhffSg7deGJh9/OOJL88Zr1fTfK2bvodSb42Zrn7k8zpfl+X5OZu+s+TfGDMnTE2jveZI8q/Lr07cBy828aJSU7uphend2eN0n3P/UnO6967J8mvd9N/ld7z8f4049xJZkxdvK1b1z1J7h+n/Hek9zirkt6jg15Ob0/0hCRbDm77KP8OO3PEsxG9vCb70p6nTXuem+TRJO8a9t/EoF4j2QPtPFNr/dtu+k+SXJLkrCRP11qf7Obfmd7DW19KsjfJ7aWUq5LsOXJlpZQ3JnlnkntLKY8luS3J2Gf33VtrPTBOOX4uyWe76T/uyvHjPpN0z7I7uPkk/62UsjW921idntce5/N0rfWxbnpLeo3woNvSeyjtjeOs/+zusztqr3X8yVHKkSR/3i3zeJLnaq2P11p/lN7typYd43PQL9rz8NvzJ5I8WGt96BjLzCqjHKBH3sPwqPc0rLXuT++pCfcleW+SvxhnsROSvFhrPW/M6+1j3n+loYzH+szY9/5DkiVJLqy1npfkufT2OJPk1THLHcjhj7D7uyQ/X0qZl+NzcBs/OmJ7P4pH5jEY2vMQ23Mp5fquzP/lOLc9o4xygP50KeXnuumrk/xNesM6y0opb+3m/0aSv+72Rk+ptX4pyR8kWdm9/69JFiRJrfXlJE+XUn4lSUrPweWO5e/Se5pF0ms4LXtvpyTZXWvdV0r5+SRnTvBztyf5UpJ7yutPXtieXl28pfv91xrKBYOiPQ+pPZdSfjvJ5Ul+reupjoxRDtBvJfndUsoTSU5NckutdW+Sa9Mbtnk8vT2uW9NrVPd3Qyp/k9f2sv4syfruJIW3pNdgfquUcvBpC2snUI7/mOTabt2/keQ/N3yXu5Os6sp8TXqNZUJqrf8jydeT/HEp5YQx8/cm+Z0kX+xOOth9lFVMWinlP5VSdqV3ssHWUsqn+7VuRpb2nOG05/Tq9E1J/k93AtZ/7eO6p7WRfBpLKWVZegfQzxl2WYDjoz0zLKPcAwWAZiPZAwWA46UHCgANBCgANBCgANBAgAJAAwEKAA0EKAA0+P8kFyVfQTrlwwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 576x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ranks, dap_samples = run_sbc(thetas, xs, posterior_)\n",
    "check_stats = check_sbc(ranks, thetas, dap_samples, 1_000)\n",
    "print(check_stats)\n",
    "f, ax = sbc_rank_plot(ranks, 1_000, plot_type=\"hist\", num_bins=30)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The histogram of ranks now shoots above the allowed (greyed) area for a uniform distributed around the extrema. We made the posterior samples too thin, so we received more extreme counts of ranks. The KS p-values vanish again and the `c2st` metric of the ranks is also larger than `.5` which underlines that our rank distribution is not uniformly distributed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We again see, **the rank distribution is capable of identifying pathologies of the estimated posterior**:\n",
    "\n",
    "+ a **centrally peaked rank distribution** shows a systematic **over-estimation of the posterior variance**  \n",
    "(we dispersed the variance of the posterior by a factor of `2`)\n",
    "+ a **U shaped rank distribution** shows a systematic **under-estimation of the posterior variance**  \n",
    "(we dispersed the variance of the posterior by a factor of `.5`)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Simulation-based calibration offers a direct handle on which pathology an estimated posterior examines. Outside of this tutorial, you may very well encounter situations with mixtures of effects (a shifted mean and over-estimated variance). Moreover, uncovering a malignant posterior is only the first step to fix your analysis. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "2193897e41726b46f35b9de052100742a934a9183b8a000ae8eb69e12e860d83"
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.11"
  },
  "name": "13_diagnosis_sbc.ipynb"
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
